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ABSTRACT 

In this paper we describe a uniform analysis of eight transits and eleven secondary eclipses of the 
extrasolar planet GJ 436b obtained in the 3.6, 4.5, and 8.0 /ini bands using the IRAC instrument on 
the Spitzer Space Telescope between UT 2007 June 29 and UT 2009 Feb 4. We find that the best-fit 
transit depths for visits in the same bandpass can vary by as much as 8% of the total (4.7cr significance) 
from one epoch to the next. Although we cannot entirely rule out residual detector effects or a time- 
varying, high-altitude cloud layer in the planet's atmosphere as the cause of these variations, we 
consider the occupation of active regions on the star in a subset of the transit observations to be the 
most likely explanation. We find that for the deepest 3.6 /xm transit the in-transit data has a higher 
standard deviation than the out-of-transit data, as would be expected if the planet occulted a star 
spot. We also compare all published transit observations for this object and find that transits observed 
in the infrared typically have smaller timing offsets than those observed in visible light. In this case 
the three deepest Spitzer transits are all measured within a period of five days, consistent with a 
single epoch of increased stellar activity. We reconcile the presence of magnetically active regions 
with the lack of significant visible or infrared flux variations from the star by proposing that the star's 
spin axis is tilted with respect to our line of sight, and t hat the planet's or bit is therefore likely to 
be misaligned. In contrast to the results reported by Beaulieu et al.l (|2011h . we find no convincing 
evidence for methane absorption in the planet's transmission spectrum. If we exclude the transits that 
we believe to be most affected by stellar activity, we find that we prefer models with enhanced CO 
and reduced methane, consistent with GJ 436b's dayside composition from I Stevenson et al.l (|2010[ ). It 
is also possible that all transits are significantly affected by this activity, in which case it may not be 
feasible to characterize the planet's transmission spectrum using broadband photometry obtained over 
multiple epochs. These observations serve to illustrate the challenges associated with transmission 
spectroscopy of planets orbiting late-type stars; we expect that other systems, such as GJ 1214, may 
display comparably variable transit depths. We compare the limb-darkening coefficients predicted by 
PHOENIX and ATLAS stellar atmosphere models, and discuss the effect that these coefficients have on 
the measured planet-star radius ratios given GJ 436b's near-grazing transit geometry. Our measured 
8 /im secondary eclipse depths are consistent with a constant value, and we place a la upper limit of 
17% on changes in the planet's dayside fiux in this band. These results are consistent with predictions 
from general circulation models for this planet, which find that the planet's dayside flux varies by a 
few percent or less in the 8 /xm band. Averaging over the eleven visits gives us an improved estimate 
of 0.0452% ± 0.0027% for the secondary eclipse depth; we also examine residuals from the echpse 
ingress and egress and place an upper limit on deviations caused by a nonuniform surface brightness 
for GJ 436b. We combine timing information from our observations with previously published data 
to produce a refined orbital ephemeris, and determine that the best-fit transit and eclipse times are 
consistent with a constant orbital period. We find that the secondary eclipse occurs at a phase of 
0.58672 ± 0.00017, corresponding to ecos(a;) = 0.13754 ± 0.00027 where e is the planet's orbital 
eccentricity and lo is the longitude of pericenter. We also present improved estimates for other system 
parameters, including the orbital inclination, a/i?*, and the planet-star radius ratio. 

Subject headings: binaries: eclipsing — stars: activity — planetary systems — techniques: photometric 
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1. INTRODUCTION 

Transiting planet systems have proven to be a powerful 
tool for studying exoplanetary atmospheres. Observa- 
tions of transiting systems have been used to detect the 
signatures of atomic and molecular absorption features 
at wa velen gths rangi ng fr oi n the UV to the i nfrared 
(e.g. J iCh arbonneau et al.l l2002t iVidal- Madi ar et al 



2003"; Swain ct al. 2008; 'Desert et al.l 120081 : iPont et al 
2008a: Linskv ct al.. .201 0). although some times the re- 



sults have proven to be controversial (e.g.. iGibson et all 
120 lit) . They have enabled studies of the dayside 
emission spectra and pre ssure-temperature p rofiles 
of close-in planet s (e.g., Charbonneau et al. 120051: 
Dem ing ct al. 2005: (Knutson ct al. 2008; Grillma ir et al.l 
2008f ). and they have inf ormed us about thei r atmo- 



spheric circulation (e.g., iKiiutson et al.l 120071 I20d9at 
Cowan. Agol. fc Charbonneaul l2007t iCrossfield et al.l 
20101) . Although we currently know of 101 transiting 



planet systems, our knowledge of these planets (includ- 
ing a majority of the studies cited above) has so far been 
dominated by studies of the brightest and closest handful 
of systems, including HD 209458b and HD 189733b. 
Planets orbiting small stars offer additional advantages, 
as they produce proportionally deeper transits and 
secondary eclipses as a result of their favorable radius 
ratios and lower stellar effective temperatures. B y this 
standard, GJ 436 (Butler ct al. 2004; Ma ness et alfcoOTt 
Gillon et al 1l2007a.b: Deming et al..,2007i nD'emorv et al.l 



2007t iTorr cs 2007) represents an ideal target, as the 



primary in this system is an early M star with a K band 
magnitude of 6.1. 

GJ 436b is currently one of the smallest known tran- 
siting p lanets, with a mass only 22 times that of the 
Earth (jTorresI l2007f l. Of the planets orbiting stars 
brighter than 9th magnitu de in K band, only GJ 1214b 
(jCharbonneau et all I2009D . which also orbits a nearby 
M dwarf, is smaller. New discoveries of low-mass tran- 
siting planets from space-based surveys such as the Ke- 
pler and CoRoT missions are unlikely to change this pic- 
ture significantly, as both include relatively few bright 
stars. GJ 436b is also one of the coolest known transit- 
ing pla nets, with a dayside ef fective temperature of only 
800 K (jStevenson et al.|[20lol ). Like GJ 1214b, GJ 436b 
has a high average density indicative of a massive rocky 
or icy core. In GJ 436b's case, models indicate that it 
must also maintain 1 — 3 Mm of its mass in the form of 
a H/He atmosphere (lAdams et a l.l '200^ iFigueira et al.l 
[20091: iRogers fc SeageHl2010l: iNettclnian n et al.ll2010D in 
order to match the observed radius. 

By measuring the wavelength-dependent transit depth 
as GJ 436b passes in front of its host star we can study 
its atmospheric composition at the day-night terminator, 
which should b e dominated by methane, water, and car- 
bon monoxide (iSpiegel e t al."2010': 'Stevenson et al."201d: 



Shabram ct al. J2011I: Madhusudhan & Scager 2011) 



Pont e t al. (2008^ observed two transits of GJ 436b with 



NICMOS grism spectrograph on the Hubble Space Tele- 
scope (HST) and placed an upper limit on the ampli- 
tude of the predicted water absorption feature between 
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1 — 2 fim. More recentlv iBeaulieu et al.l ()201lD reported 
the detection of strong methane absorption in the 3.6, 
4.5, and 8.0 /im Spitzer bands. 

We can compare these results to observations of the 
planet's dayside emission spectrum, obtained by measur- 
ing the depth of the s econdary eclipse when the planet 
passes behind the star. iStevenson et al.l ()2010f ) measured 
secondary eclipse depths for GJ 436b in the 3.6, 4.5, 5.8, 
8.0, 16, and 24 /im Spitzer bands, from which they con- 
cluded that the planet's dayside atmosphere contained 
significantly less methane and more CO than the equi- 
librium chemistry predictions. In this work we present 
an analysis of eight transits and eleven secondary eclipses 
of GJ 436b observed with Spitzer, including an indepen- 
dent analysis of the transit data described in Beaulieu 
et al., and discuss the corresponding implications for GJ 
436b's atmospheric composition. 

Unlike most close-in planets, which typically have cir- 
cular orbits, GJ 43 6b has an orbital eccentricity of 
approximatelv 0.15 ([Maness et al.' I 2007t iDem ing et al.l 
120071: iDemorv et al.ll2007l: [Madhusud han fc Wi nn 2009) . 
Atmospheric circulation models for eccentric Jovian 
planets suggest that they may exhibit significant 
temperature variations fro n i one orbit t o the next 
(iLangton fc LaughlinI [2008t liro fc Demiil^ [20Tol) . al- 
though [Lewis et al.l ([20101 ) find little evidence for sig- 
nificant temporal variability in general circulation mod- 
els for GJ 436b. The extensive nature of our data set, 
which includes eleven secondary eclipse observations in 
the same bandpass obtained between 2007 — 2009, allows 
us to test the predictions of these models by searching 
for changes in the planet's 8 fim. dayside emission on 
timescales ranging from weeks t o years. 

It has also been suggested ([Ribas et al.l [lOOS") that 
GJ 436b's orbital parameters are changing in time, per- 
haps through perturbations by an unseen second planet 
in the system. Such a planet could serve to maintain 
GJ 436b's nonzero eccentricity despite ongoing orbital 
circularization, and would not necessarily produce transit 
timing variations large enough to be detected by earlier, 
ground-based studi es (iBatvg in et al. 20 091). Although 



120081: ICoughlin et all 


20081: IMadhusudhan & Winnll2009l: 


[Caceres et al.l [20091 


Shoorer et al.l [20091: [Ballard et al.l 



varying orbital parameters or a second transiting object 
in the system, Spitzer^ s unparalleled sensitivity and sta- 
bility allow us to extend the current baseline by nine 
months with new, high-precision transit observations. 

2. OBSERVATIONS 

We analyze nineteen separate observations of GJ 436, 
including two 3.6 fim transits, two 4.5 ^m transits, four 
8 /im transits, and eleven 8 /im secondary eclipses, as 
listed in Table [T] All observations were obtai ned between 
2007 and 2009 using the IRAC instrume nt ([Fazio et al.l 
|2004[) on the Spitzer Space Telescope (Wer ner et al.ll2004D 
in subarray mode. Some of these data were previ- 
ously published by other groups, including a transit 
and secondary eclipse observed on UT 2007 Jun 29/30 
(IDeming et al.l l2007t iGillon et all l2007al: IDemor v et aT 



2007 ) and six transits observed in 2009 (fBcaulicu et al 



20111) . Because the two shorter wavelength IRAC chan- 



nels (3.6 and 4.5 /xm) use InSb detectors and the two 
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Table 1 

Spitzer Observations of GJ 436b 



UT Date 


Event 


A (^m) 


Duration (hr) 


tint (s) 


'^exposures 


Bkd (MJy Sr-l)'' 


Flux (MJy Sr-l)'*''^ 




U i zUUY Jun zy 


Transit 


8.0 


3.4 


0.4 


oo A on 

28,480 


530.8 


9,149.0 


0.50070 


U i zUU7 Jun oU 


Eclipse 


8.0 


5.9 


0.4 


49,920 


544.2 


9,148.2 


0.49870 


U i 2008 Jun 11 


Eclipse 


8.0 


3.4 


0.4 


28,800 


226.0 


9,156.6 


0.49770 


U i 2UUo Jun lo 


Eclipse 


8.0 


3.4 


0.4 


28,800 


255.0 


9,161.3 


0.502% 


U 1 iOOs Jun Id 


Eclipse 


8.0 


3.4 


0.4 


28,800 


296.5 


9,154.5 


0.494% 


U 1 200o Jun 19 


Eclipse 


8.0 


3.4 


0.4 


28,800 


306.4 


9,151.7 


0.493/0 


UT 2008 Jul 12^^ 


Eclipse 


8.0 


70 


0.4 


588,480 


669.2 


9,159.9 


0.506% 


UT 2008 Jul 14^^ 


Transit 


8.0 


70 


0.4 


588,480 


695.6 


9,160.4 


0.509% 


UT 2008 Jul IS'^ 


Eclipse 


8.0 


70 


0.4 


588,480 


714.4 


9,158.1 


0.515% 


UT 2009 Jan 9 


Transit 


3.6 


4.3 


0.1 


117,056 


37.7 


36,164.3 


0.387% 


UT 2009 Jan 17 


Transit 


4.5 


4.3 


0.1 


117,056 


61.6 


24,382.6 


0.561% 


UT 2009 Jan 25 


Transit 


8.0 


4.3 


0.4 


35,904 


474.5 


9,151.9 


0.502% 


UT 2009 Jan 27 


Eclipse 


8.0 


3.4 


0.4 


28,800 


455.1 


9,164.5 


0.495% 


UT 2009 Jan 28 


Transit 


3.6 


4.3 


0.1 


117,056 


82.5 


36,744.5 


0.389% 


UT 2009 Jan 29 


Eclipse 


8.0 


3.4 


0.4 


28,800 


411.8 


9,161.7 


0.496% 


UT 2009 Jan 30 


Transit 


4.5 


4.3 


0.1 


117,056 


86.3 


24,177.0 


0.567% 


UT 2009 Feb 1 


Eclipse 


8.0 


3.4 


0.4 


28,800 


395.9 


9,163.9 


0.499% 


UT 2009 Feb 2 


Transit 


8.0 


4.3 


0.4 


35,904 


393.6 


9,143.8 


0.501% 


UT 2009 Feb 4 


Eclipse 


8.0 


3.4 


0.4 


28,800 


376.6 


9,154.5 


0.499% 



^ Average sky backgrounds and stellar fluxes estimated for a 5 pixel aperture. 

^ In order to minimize the effects of the detector ramp in the 8.0 pm observations wc estimate the out-of-transit flux using data after the end 
of the eclipse event where the ramp is generally smallest; for consistency wc use the same region to estimate the fluxes at 3.6 and 4.5 /im. We 
use a 5.0 pixel aperture for the photometry at [3.6,4.5,8.0] ^m and apply the appropriate aperture correction of [1.049,1.050,1.068] from Table 
4.7 of the IRAC Instrument Handbook to determine the total flux from the star in each observation. 
^ Standard deviation of residuals after dividing out best-flt corrections for instrument effects and transit light curves. 

These events were observed as part of a single, continuous phase curve observation with a duration of 70 hours spanning two secondary 
eclipses and one transit. 



longer wavelength channels (5.8 and 8.0 /im) use Si:As 
detectors, each of which display different detector effects, 
we describe our analysis for each type separately below. 

We calculate the BJD_UTC values at mid-exposure for 
each image using the DATE_OBS keyword in the im- 
age headers and the position of Spitzer^ which is in an 
earth-trailing orbit, as determined using the JPL Hori- 
zons ephemeris. Each set of 64 images obtained in subar- 
ray mode comes as a single FITS file with a time stamp 
corresponding to the start of the first image; we calculate 
the time stamps for individual images assuming uniform 
spacing and using the difference between the AINTBEG 
and ATIIvIEEND headers, which record the start and 
end of the 64-image series . We t hen use the routines de- 
scribed in lEastman et all (|201Clf ) to convert from Spitzer 
JD to BJD_UTC. Eastman et al. further advocate a con- 
version from UTC to TT timing standards, which pro- 
vide a more consistent treatment of leap seconds. We 
note that for the dates spanned by these observations the 
conversion from BJD_UTC to BJD_TT simply requires 
the addition of [65.184,65.184,66.184] s for data obtained 
in [2007,2008,2009], and we leave the dates listed in Table 
[5] in BJD_UTC for consistency with other studies. 

2.1. 3.6 and 4-5 fJ-m Photometry 

GJ 436 has a K band magnitude of 6.07, and as a 
result we elect to use short 0.1 s exposures at 3.6 and 
4.5 /im in order to ensure that we remain well below sat- 
uration. Subarray images have dimensions of 32 x 32 
pixels, making it challenging to estimate the sky back- 
ground independent of contamination from the wings of 
the star's point spread function. We choose to exclude 
pixels within a radius of 12 pixels of the star's position, 
as well as the 14th- 17th rows, which contain a horizon- 
tal diffraction spike that extends close to the edges of 



the array. We also exlcude the top (32nd) row of pix- 
els, which have values that are consistently lower than 
those for the rest of the array. We then iteratively trim 
3cr outliers from the remaining subset of approximately 
six hundred pixels, create a histogram of the remaining 
values, and then fit a Gaussian to this histogram to de- 
termine the sky background for each image. We find that 
the background is 0.1 - 0.2% and 0.3 - 0.4% of the total 
fiux in a 5 pixel aperture for the 3.6 and 4.5 /im arrays, 
respectively. 

We correct for transient hot pixels by taking a 10- 
pixel running median of the fiuxes at a given pixel posi- 
tion within each set of 64 images and replacing outliers 
greater than Aa with the median value. We found that 
using a wider median filter or tighter upper limit for dis- 
criminating outliers increased the scatter in the final time 
series while failing to significantly reduce the number of 
images that are ultimately discarded. We find that ap- 
proximately 0.4 — 0.8% of our images have one or more 
pixels fiagged as outliers using this filter. 

Several recent papers have investigated optimal meth- 
ods for estimating the position of the star on the array 
for Spitzer photo metry, with the most e xten sive discus- 
sions appearing in lStevenson et all (|2010[ ) and Agol et al.l 
(pOlOi) . These papers conclude that flux-w eight ed cen- 
troiding (e.g., IKnutson et al.l [20081 : ICharbonneau et al. 
|200^ and parabola-fitting routines (e.g.. iDeming et ah 
2006, 2007) tend to produce less than optimal results, 
while Gaussian fits and least asymmetry methods ap- 
pear to have fewer systematic biases and a lower overall 
scatter. We confirm that for all three wavelengths we ob- 
tain better results (defined as a lower scatter in the final 
trimmed light curve after correcting for detector effects) 
with Gaussian fits than with flux-weighted centroiding, 
with a total reduction of 2 — 7% in the standard deviation 



4 



Knutson et al. 




Figure 1. Raw photometry for the eight observed transits of GJ 436b, arranged in chronological order and with best-fit detector functions 
overplotted (solid red lines). Data has been binned in either 0.9 minute (3.6, 4.5 ^m) or 1.5 minute (8.0 fim) bins. 



of the final time series binned in sets of 64 images. We 
obtain the best results in both the 3.6 and 4.5 /im bands 
when we first subtract the best-fit background flux from 
each image, correct bad pixels as described above, and 
then fit a two-dimensional Gaussian function to a cir- 
cular region with a radius of 4 pixels centered on the 
position of the star. Using smaller or larger fitting re- 
gions does not significantly alter the time series but does 
result in a slightly higher scatter in the normalized light 
curve. Although error arrays are available as part of the 
standard Spitzer pipeline, we find that in this case we 
obtain better results using uniform error weighting for 
individual pixels. We use a radially symmetric Gaussian 
function and run our position estimation routines twice, 
once where the width is allowed to vary freely in the fits 
and a second time where we fix the width to the me- 
dian value over the time series. Reducing the degrees of 
freedom by fixing the width produces fits that converge 
more consistently, with a corresponding improvement in 
the standard deviati on of the normali z ed tim e series and 
fewer large outliers. iStevenson et al.l (|201Cll ) report that 
they obtain better position estimates when fitting Gaus- 
sians to images that have been interpolated to 5 x higher 
resolution, but we find that using interpolated images for 
our position fits resulted in a slight increase in the scatter 
in our final light curves. 
We perform aperture photometry on our images using 



the position estimates derived from our Gaussian fits; 
we expect that aperture photometry will produce the 
optimal results in light of the low background flux at 
these wavelengths. We use apertures with radii ranging 
between 2.5-7.0 pixels in half pixel steps. We find that 
apertures smaller than 3.5 pixels show excess noise, likely 
connected to position-dependent flux losses, while aper- 
tures larger than 5 pixels are more likely to include tran- 
sient hot pixels and higher background levels, resulting 
in a higher root-mean-square variance in the flnal light 
curve. We use a 3.5 pixel aperture for our final analy- 
sis, but we find consistent results for apertures between 
3.5 — 5.0 pixels. We trim outliers from our final time 
series using a 50 point running median, where we dis- 
card outliers greater than 3a, approximately 2% of the 
points in a typical light curve. We find that we trim 
fewer points when we use flux-weighted centroiding for 
our position estimates (typically 0.6%), but the uncer- 
tainties in our best-flt transit parameters are still larger 
than with the Gaussian fits due to the increased scatter 
in the final trimmed time series. We also trim the first 15 
minutes in all observations except for the 4.5 /xm transit 
on UT 2009 Jan 30, where we trim the first hour of data. 
Images taken at the start of a new observation tend to 
have larger pointing offsets, most likely due to the set- 
tling time of the telescope at the new pointing; we find 
that discarding these early data improves the quality of 
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the fit to the subsequent points. For all visits other than 
the transit on UT 2009 Jan 30, we find that we achieve 
consistent results when we trim either the first 15, 30, 
or 60 minutes of data, and we therefore choose to trim 
the minimal 15-minute interval. For the UT 2009 Jan 
30 observation we find that the data display an addi- 
tional time dependence that is not well-described by the 
standard linear function of time in Eq. [Jl but is instead 
better-described by a linear function of ln((ii). This may 
be due to the fact that the star falls near the edge of a 
pixel in these observations, which could introduce addi- 
tional time-dependent effects. Rather than changing the 
functional form used to fit these data, we instead opt to 
trim the first hour of observations, which removes the 
most steeply- varying part of the time series and leaves a 
trend that is well-described by the same linear function 
of time used in the other transit fits. We find that we 
obtain the same transit parameters for this visit when 
we either trim the first 15 minutes of data and fit with 
a linear function of In(dt) or trim the first hour of data 
and fit with a linear function of time, so this choice does 
not affect our final conclusions. 

Fluxes measured at these two wavelengths show a 
strong correlation with the changing position of the 
star on the array, at a level comparable to the 
depth of the secondary eclipse. This effect is due 
to a well-documented intra-pixel sensitivity variation 
(e.g., [R each et al. 2005; Charbonneau ct al. 2005, 200| 
iMorales-Calderon et al.l I2006t IKnutson et al., 2008), in 
which the sensitivity of an individual pixel differs by 
several percent between the center and the edge. The 
3.6 /im array typically exhibits larger sensitivity varia- 
tions than the 4.5 /im array, as demonstrated by the UT 
2009 Jan 9 and 17 transits. The UT 2009 Jan 30 transit 
falls very near the edge of a pixel in the 4.5 /im subar- 
ray, and thus displays a sensitivity variation comparable 
to that of the more centrally-located 3.6 /xm transit on 
UT 2009 Jan 28. We correct for these sensitivity varia- 
tions by fitting a quadratic function of x and y position 
simultaneous with the transit light curve: 

/= fQ*{ai+a2{x-xo) + a'i{x-xaf 

+ai{y - Vq) + a<i{y - yaf + aat) (1) 

where /o is the original flux from the star, / is the 
measured flux, x and y denote the location of the star 
on the array, xq and yo are the median x and y posi- 
tions, t is the time from the predicted eclipse center, and 
oi — og are free parameters in the fit. In both band- 
passes we find that quadratic terms in both x and y are 
necessary to achieve a satisfactory fit to the observed 
variations, although the value for the fits is not im- 
proved by the addition of an xy term, or higher-order 
terms in x and y. We find that the fits are also im- 
proved by the addition of a linear term in time, consistent 
with previous observations at these wavelengths (e.g., 
Knutson et"alll2009bt iTodorov et al.ll201rt 'Fressin et alJ 
2010; O'Donovan et al.ll2010l : IPa ning ct al. 2010). 

2.2. 8.0 fim Photometry 

We follow the same methods described in ij2.1l to esti- 
mate the sky background in the 8.0 /im images, except 
in this case we include pixels at distances of more than 



ten pixels from the position of the star in our estimate 
instead of the previous twelve-pixel radius. The back- 
ground in these images ranges between 2.6 — 7.7% of the 
total flux in a five pixel aperture, and we find that in- 
cluding pixels between 10 — 12 pixels from the star's posi- 
tion improves the accuracy of our background estimates 
without adding significan t contamination fr om the star's 
point spread function. In lAgol et all ()2010fl we find that 
using a slightly larger 4.5 pixel aperture instead of 3.5 
pixels minimizes correlated noise in 8 /im Spitzer obser- 
vations (albeit at the cost of slightly higher Gaussian 
noise), and we therefore elect to use a 4.5 pixel aperture 
for our 8 /xm data. Our choice of aperture has a negli- 
gible effect on the best-fit eclipse depths and times, as 
we find consistent results for apertures between 3.5 — 5.0 
pixels. 

Spitzer fiuxes for stars observed using the IRAC 
8.0 iim array, the IRS 16 /im array, and the MIPS 
24 /jm array do not appear to have a significant posi- 
tion dependence, but do display a ramp-like behavior 
where higher-illumination pixels converge to a constant 
value within the first hour of observations while lower- 
illumination pixels show a continually increasing linear 
trend on the time scales of interest here. This effect is 
believed to be due to charge-trapping in the array, and 
is discussed in de tail in lKnutson et al.l (|2007l [2009c) and 
lAeol et al.l(l20T0l) . among others. We mitigate this effect 
in our data by staring either at a bright star (HD 107158 
in the case of the 8 fim secondary eclipse observations 
between UT 2008 Jun 11 and Jun 19) or an HII region 
with bright diffuse emission at 8 /im (LBN 543 for the 
8 /im transit observations, and Gill. 612-1-0. 374 for the 
8 /im secondary eclipse observations between UT 2009 
Jan 27 and Feb 4) for approximately 30 minutes prior 
to the start of our observations. The 2007 observations 
were obtained prior to the de velopment of this pre flash 
technique, but as discussed in iDeming et al.l ()2007D the 
transit observation happened to follow an observation of 
another bright object and thus was effectively prefiashed 
in the same manner as the 2008 and 2009 data. The sec- 
ondary eclipse observed in 2007 was not preflashed, and 
thus displays a much steeper ramp than the other obser- 
vations. We examine the distribution of ramp slopes in 
Fig. Uland find no correlation between the relative offsets 
in the positions of GJ 436 and the prefiash star and the 
slope of the subsequent ramp; the prefiash star is offset 
by [0.4, 0.2, 0.3, 0.1] pixels in the UT 2008 Jun 11, 13, 
16, and 19 observations, respectively, but the shallowest 
ramp occurs in the Jun 13 observation while the small- 
est offset occurs in the Jun 19 observation. We speculate 
that the UT 2008 Jun 13 observation may have been ef- 
fectively preflashed by the preceding science observations 
in the same way as the UT 2007 Jun 29 transit observa- 
tion. We find that all forms of prefiash reduce the slope of 
the subsequent ramp as compared to the non-prefiashed 
secondary eclipse on UT 2007 June 30, but the HII re- 
gions consistently produce a larger reduction in the ramp 
slope than preflashes using a bright star. 

We can describe the ramps in our 8 fim science data 
with the following functional form; 

/ = ci (1 - C2 exp {-St/cs) - C4 exp(-(5t/ C5)) (2) 
where / is the measured flux, 6t is the elapsed time from 
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Figure 2. Photometry for the eight observed transits of GJ 436b after the best-fit corrections for instrument effects are removed, arranged 
in chronological order. Data has been binned in either 2.7 minute (3.6, 4.5 ^tm) or 4.3 minute (8.0 fim) bins. Best-fit transit curves are 
overplotted in red, and the residuals from each fit are shown in the lower panel. In this plot we have assumed a constant ephemeris for 
the planet rather than using the best-fit transit times. Note that although the out-of-transit residuals for the second 3.6 fim observation 
on UT 2009 Jan 28 appear to be relatively Gaussian, there are additional variations during the transit that are not well-accounted for 
by the best-fit transit light curve. These variations are likely due to occultations of spots or faculae by the planet. The residuals for the 
4.5 fim transit observed on UT 2009 Jan 30 display excess correlated noise both in and out of transit, most likely due to an imperfect 
correction for the sharp flux variations caused by the star's location at the edge of a pixel. 



the start of the observations, and ci — C5 are free pa- 
rameters in the fit. Previous studie s have elected to us e 
either a single exponential (e.g.. iHarrington et al.ll2007[ ). 
a linear -I- log function of St (e.g.. lDerning et al.ll2007D. or 
a quadratic function in lo^j St) (e.g., Charbo nneau et al] 
[2OO8; Knutson ct al. 2008; Desert et a l. 2009). However, 
in lAgol et al.l ()2010i ) we find that the functional forms in- 
volving log((5t) produce eclipse depths that are correlated 
with the slope of the observed ramp function, while the 
single exponential does not provide a good fit to data 
with a steep ramp. We conclude that a double expo- 



nential function has enough degrees of freedom to fit a 
range of ramp profiles, while still avoiding correlations 
between the measured eclipse depths and the slope of 
the detector ramp. Although we require a double expo- 
nential function in order to fit the steeper, non-preflashed 
2007 secondary eclipse observation in this study, we ob- 
tain comparable results with a single exponential term 
for our prefiashed 8 data. We therefore elect to use 
this simpler single exponential in our subsequent analysis 
for all 8 /im visits except the 2007 secondary eclipse. 
For our fits to phase curve data obtained on UT 2008 




Figure 3. Standard deviation of residuals vs. bin size for the eight transits observed with Spitzer, arranged in chronological order. 
Observations wrere taken at 8.0, 8.0, 3.6, 4.5 fim (top row, left to right), 8.0, 3.6, 4.5, and 8.0 /im (bottom row, left to right), respectively. 
The red curve shows the predicted root-n scaling expected for Gaussian noise. 



Table 2 

Individual Best-Fit Transit Parameters 



UT Date 



Depth 



Transit Center (BJD) O-C (s)'' 



UT 2007 Jun 29 8.0 0.08322 ± 0.00052 0.6926% ± 0.0087% 

UT 2008 Jul 14 8.0 0.08247 ± 0.00061 0.6801% ± 0.0101% 

UT 2009 Jan 9 3.6 0.08182 ± 0.00037 0.6694% ± 0.0061% 

UT 2009 Jan 17 4.5 0.08286 ± 0.00047 0.6865% ± 0.0078% 

UT 2009 Jan 25 8.0 0.08224 ± 0.00051 0.6763% ± 0.0084% 

UT 2009 Jan 28 3.6 0.08495 ± 0.00056 0.7216% ± 0.0095% 

UT 2009 Jan 30 4.5 0.08502 ± 0.00057 0.7227% ± 0.0097% 

UT 2009 Feb 2 8.0 0.08424 ± 0.00049 0.7096% ± 0.0083% 



2454280.78193 ± 0.00012 
2454661.50314 ± 0.00017 
2454841.28821 ± 0.00008 
2454849.21985 ± 0.00012 
2454857.15155 ±0.00012 
2454859.79504 ± 0.00012 
2454862.43970 ± 0.00011 
2454865.08345 ± 0.00012 



12.5 ± 10.2 
6.2 ±14.4 
6.9 ± 6.5 

2.1 ± 10.5 

3.2 ±10.1 
-31.9 ± 10.3 

33.7 ±9.6 

20.8 ± 10.6 



^ Observed minus calculated transit times. Predictions use the best-fit ephemeris of = 2454865.083208 ± 0.000042 BJD 
and P = 2.6438979 ± 0.0000003 days from Table[5] 



July 12-15, we select a four hour subset of data centered 
on the position of the transit or eclipse and use that 
in our fits. The first eclipse takes place at the start of 
the observations, which exhibit a residual ramp, and we 
therefore fit this light curve with the same single expo- 
nential as our other data. We use a linear function of time 
to fit the out-of-eclipse trends in the transit, which occurs 
in the middle of the observations, as well as the secondary 
eclipse towards the end of the observations. We find that 
the scatter in the central region of the time series near the 
transit, when the star is closest to the edge of the pixel. 



is higher than for either sec ondary eclipse or for th e other 
8 /im transit observations. ([Stevenson et al.ll2010t ) found 
that the fluxes measured with the 5.8 fiva. Spitzer array 
sometimes display a weak dependence on the position of 
the star, which may be due to either flat-fielding errors or 
intrapixel sensitivity variations similar to those observed 
in the 3.6 and 4.5 /im arrays, although no such effect 
has been definitively detected in the 8 /im array to date. 
We test for the presence of position-dependent flux vari- 
ations in our data by adding linear functions of x and y 
position to each of our 8 /Ltm transit fits, and find that the 
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Figure 4. Raw photometry for eleven 8 /im secondary eclipses of GJ 436b, arranged in chronological order. Data has been binned in 2.2 
minute bins, and the best-fit corrections for detector effects in each visit are overplotted in red. 



value of the resulting fits is effectively unchanged in 
all cases except for the UT 2008 July 14 transit, where it 
decreased from 37,186.6 to 37,177.7 for 33,636 points and 
six degrees of freedom. Using the Bayesian Information 
Criterion described in Stevenson et al. (2010), we con- 
clude that this reduction in is not significant, and we 
exclude these position-depedent terms in our subsequent 
analysis of the 8 /im data. 

2.3. Transit and Eclipse Fits 

We carry out simultaneous fits to determine the best- 
fit transit functions and detector corrections using a 
non-linear Levenb erg-Marquardt minimization routine 
(|Markwardtl 120091 ) . We calcula te our eclipse curve us- 
ing the equations from jMandcl fc Agoll ([2002:) assuming 
a longitude of pericenter equal to 334° ± 10° (update 
based on complete set of published and unpublished ra- 
dial velocity data, A. Howard, personal communication, 
2010). The orbital eccentricity determined from the up- 
dated radial velocity data is 0.145 ±0.017, but we choose 
to set the orbital eccentricity equal to 0.152 in our fits, 
which we calculate using the above longitude of pericen- 
ter a nd the published value o f e*cos{uj) = 0.1368±0.0004 
from iStevenson et al.l (|2010[ ). We find that the uncer- 
tainty in the calculated eccentricity is dominated by the 
uncertainty in w, but this has a minimal impact on our 
transit fits. Our best-fit parameters change by less than 



la for eccentricity values between 0.142 and 0.169, cor- 
responding to ±10° in w, where our best-fit inclination is 
most sensitive to the assumed eccentricity (0.9ct change), 
a/Ri, is somewhat sensitive (0.5(t change), and the best- 
fit radius ratios and transit times for individual fits are 
minimally sensitive (< 0.3a change). Our nominal values 
for the eccentricity and longitude of pericenter result in 
a transit length of 60.9 minutes, 0.5 minutes longer than 
the zero eccentricity case. Using the same parameters for 
the secondary eclipse, which occurs shortly before peri- 
astron passage, produces a length of 62.6 minutes. 

We fit the eight transits simultaneously and assume 
that the inclination and the ratio of the orbital semi- 
major axis to the stellar radius o/i?+ are the same for 
all transits, but allow the planet-star radius ratio Rp/Ri, 
and transit times to vary individually. Figure[T]shows the 
final binned data from these fits with the best fit normal- 
izations for the detector effects and transit light curves in 
each channel overplotted, and Figure [2] shows the binned 
data once these trends are removed, with best-fit tran- 
sit curves overplotted. Best-fit parameters are given in 
Tables [Hand El 

2.3.1. A Comparison of ATLAS and PHOENIX 
Limb-Darkening Models 

We derive limb-darkening coefficients for the star us- 
ing a Kurucz ATLAS stellar atmosphere model with 
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Figure 5. Photometry for eleven 8 secondary eclipses of GJ 436b, arranged in clironological order. Data has been binned in 6.4 
minute bins, and the best-fit eclipse curve for each visit is overplotted in red. The residuals for each visit are shown in the panels below 
the eclipses. 



r^ff = 3500 K, l ogfq) = 5.0, and [Fe/H]= (IKuruczl 
119791 11994 I2005D . where we take the flux- weighted av- 
erage of the intensity profile in each IRAC band and 
then fit this profile with four-parameter nonlinear limb- 
darkening coefficients (Claret 2000). We also derive limb- 
darkening coefficients f or a PHOENIX atmosphere model 
(|Hauschildt et allll999f) with the same parameters, and 
list both sets of coefficients in Table ID We trim the max- 
imum stellar radius in the PHOENIX models, which is set 
to an optical depth of 10~^, to match the level of the 
r = 1 surface in each Spitzer band. We estimate the lo- 
cation of this surface by determining when the intensity 
relative to that at the center of the star first drops below 
e~^, and find that the new stellar radius is 0.09 — 0.10% 
smaller than the old t = 10^^ value. We find that we can 



achieve satisfactory four-parameter nonlinear fits to the 
PHOENIX intensity profiles only when we exclude points 
where /i < 0.025, whereas the ATLAS models are well- 
described by fits including this region. 

As fllustrated in Fig. [71 the PHOENIX model predicts 
stronger limb darkening in all bands as compared to 
the ATLAS model, with the largest differences in the 
3.6 (Um band. When we compare our best-fit transit 
parameters using either the ATLAS or PHOENIX limb- 
darkening coefficients, we find that the best-fit planet- 
star radius ratios are 0.8 — 1.2ct (0.5 — 0.6%) deeper in 
the 3.6 /im band, 0.06 - 0.07cr (0.04 - 0.05%) smaller in 
the 4.5 ^m band, and 0.3 - 0.4a (0.2 - 0.3%) larger in the 
8.0 jim band for the PHOENIX models. The best-fit values 
for the inclination and a/i?* increase by 1.0a (0.6%) and 
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Figure 6. Standard deviation of residuals vs. bin size for the eleven secondary eclipses observed with Spitzer, arranged in chronological 
order. The red curve shows the predicted root-n scaling expected for Gaussian noise, and the lack of any excess noise for large bin sizes 
suggests that those light curves should be well-described by the standard MCMC error analysis. 



Table 3 

Individual Best-Fit 8 fim Secondary Eclipse Parameters 



UT Date Depth (%) 



UT 


2007 


Jun 30 





0553 







0083 


UT 


2008 


Jun 11 





0506 


± 





0110 


UT 


2008 


Jun 13 





0395 


± 





0097 


UT 


2008 


Jun 16 





0497 


± 





0087 


UT 


2008 


Jun 19 





0368 


± 





0089 


UT 


2008 


Jul 12 





0523 


± 





0090 


UT 


2008 


Jul 15 





0422 


± 





0078 


UT 


2009 


Jan 27 





0386 


± 





0087 


UT 


2009 


Jan 29 





0491 


± 





0088 


UT 


2009 


Feb 1 





0398 


± 





0086 


UT 


2009 


Feb 4 





0441 


± 





0087 



Eclipse Center (BJD) O-C (min)=' 



2454282 


3329 


± 





0016 


-0 


2 


± 


2 


3 


2454628 


6850 


± 





0017 


2 





± 


2 


4 


2454631 


3281 


± 





0021 





9 


± 


3 





2454633 


9716 


± 





0013 





3 


± 


1 


9 


2454636 


6162 


± 





0021 


1 


2 


± 


3 





2454660 


4112 


± 





0019 


1 


2 


± 


2 


8 


2454663 


0537 


± 





0040 


-0 


9 


± 


5 


8 


2454858 


7047 


± 





0026 


2 


8 


± 


3 


8 


2454861 


3460 


± 





0015 


-1 





± 


2 


2 


2454863 


9889 


± 





0017 


-2 


4 


± 


2 


4 


2454866 


6355 


± 





0023 


1 


4 


± 


3 


3 



^ Observed minus calculated transit times. Predictions use the best-fit ephemcris 
of Te = 2454865.083208 ± 0.000042 B,ID and P = 2.6438979 ± 0.0000003 days, 
and an orbital phase of 0.58685 ± 0.00017 from Tabled 
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Table 4 

Four-Parameter Nonlinear Limb-Darkening 
Coefficients'' 



Model 


Band (^tm) 


Cl 


C2 


C3 


C4 


ATLAS 


3.6 


1.122 


-1.852 


1.675 


-0.582 


ATLAS 


4.5 


0.749 


-0.917 


0.718 


-0.230 


ATLAS 


5.8 


0.815 


-1.147 


0.947 


-0.310 


ATLAS 


8.0 


0.770 


-1.141 


0.942 


-0.304 


PHDENIX 


3.6 


1.284 


-1.751 


1.433 


-0.470 


PHDENIX 


4.5 


1.203 


-1.796 


1.512 


-0.500 


PHOENIX 


5.8 


0.918 


-1.264 


1.064 


-0.358 


PHDENIX 


8.0 


0.619 


-0.762 


0.645 


-0.220 



^ Both models assume T^ff ^ 3500 K and [Fe/H]^0. The 
ATLAS model uses log(p) — 5.0 and the PHOENIX model uses 
\og{g) — 4.76 for better consistency with the radius and 
luminosity in iTorre^ | |2007l) , but emipirical tests show the 
assumed surface gravity has a negligible effect on the result- 
ing limb-darkening profiles. For a definition of this limb- 
darkening law, see 

0.9cr (0.04%), respectively, for the PHDENIX model fits; 
this is a product of the stronger limb-darkening profile, 
as GJ 436b's relatively high impact parameter creates 
a partial degeneracy between the limb-darkening profile 
and the other transit parameters. 

We examine the relative importance of the assumed 
stellar parameters by comparing two PHDENIX models 
with effective temperatures of 3400 K and 3600 K. We 
find that for this 200 K range in effective tempera- 
ture, the best-fit planet-star radius ratios change by 
0.11 - 0.16cr at 3.6 fim, 0.07 - 0.09cr at 4.5 iim, and 
0.10 — 0.12cr at 8.0 fim. The changes in the best-fit val- 
ues for the inclination and a/i?+ were similarly small, 
0.04(7 and 0.4cr, respectively. We therefore conclude that 
changes in the stellar effective temperature of less than 
200 K are negligible for the purposes of our transit fits. 
We also compute PHDENIX model intensity profiles for 
0.0 <[Fe/H]< +0.3, but we find that the differences be- 
tween models are much smaller than for our 200 K change 
in the effective temperature. 

As there are currently few observational constraints 
on liml3-darkenin!g profiles for main-sequence stars (e.g., 
IClaretl[2008L[2009ri . and even fewer constraints for M stars 
in the mid-infrared, we also consider simultaneous tran- 
sit fits in which we allow quadratic limb-darkening coeffi- 
cients in each band to vary as free parameters. As a result 
of the planet's high impact parameter, our observations 
do not directly constrain the limb-darkened intensity for 
values of 6 < 50°, corresponding to n > 0.64, as the 
planet does not cross this region on the star. However, 
we can infer the limb-darkening profile in this region if 
we assume a simple quadratic limb-darkening law. We 
require the intensity profile computed from these coef- 
ficients to be always less than or equal to one (i.e., no 
limb brightening), and we require the relative intensity 
at the edge of the star to be greater than or equa l to th e 
equivalent K band limb-darkening from iClaretl ()2000D . 
The dotted lines in Fig. [7] show the resulting best-fit 
limb-darkening profiles in each band; these profiles show 
less contrast than either model, but the ATLAS models 
appear to provide the closest match. 

This agreement is reflected in the values for the 
simultaneous transit fits; the total for the best-fit 
quadratic coefficients is 536,729.25, for the 3500 K ATLAS 
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Figure 7. A comparison of the model limb-darkening as a function 
of fj, = cos{9), where 9 ranges from 0° at the center of the star to 
90° at the edge. We show limb-darkening from four-parameter 
nonlinear limb-darkening coefficients obtained using ATLAS (solid 
lines) or PHOENIX (dashed lines) stellar atmosphere models, as well 
as the best-fit quadratic limb-darkening obtained by a fit to our 
data (dotted lines). The color indicates the bandpass, including 
3.6 /im (blue), 4.5 /im (green), and 8.0 fira (red) Spitzer bands. 
The grey shaded region indicates values of /i for which we have no 
direct observational constraints, as the planet does not cross this 
part of the star during its transit. 

limb-darkening coefficients it is 536,733.98, and for the 
[3400, 3500, 3600] K PHOENIX models it is [536,740.69, 
536,739.75, 536,738.81], for 536,798 points and either 53 
(with fixed limb- darkening) or 59 (with freely varying 
quadratic limb-darkening coefficients) free parameters. 
We use the ATLAS limb-darkening coefficients in our sub- 
sequent analysis, as they produce a marginally better 
agreement with the best fit profiles than the PHOENIX 
models. Although the value for the best-fit quadratic 
limb-darkening coefficients is formally smaller than that 
of either model, this fit also contains six additional de- 
grees of freedom, making the difference negligible. 

As an additional test, we also repeat our fits with 
the limb-darkening coefficients fixed to zero in all bands. 
This produces planet-star radius ratios that are 1.6 — 2.4cr 
(1.1%) smaller in the 3.6 /.tm band, 1.7 - 2.0cr (1.1%) 
smaller in the 4.5 fim band, and 0.9 — l.lfi (0.6 — 0.7%) 
smaller in the 8.0 /im band. The best-fit inclination and 
a/Ri, are 2.7a (1.5%)and 2.0f7 (0.1%) smaller, respec- 
tively. The value for this fit is 536,738.04, equivalent 
to the PHDENIX model fits and marginally worse than the 
ATLAS models or the fitted limb-darkening coefficients. 
This fit confirms the pattern suggested earlier, namely 
that stronger limb darkening leads to larger planet-star 
radius ratios and larger values for the inclination and 
a/Ri,. If we consider the constraints imposed by the 
transit fits, stronger limb darkening means that for a 
grazing transit the planet must occult a relatively larger 
fraction of the star in order to produce the same appar- 
ent transit depth. This effect will be even larger in visi- 
ble light, and we conclude that accurate limb-darkening 
coefficients are essential when calculating the planet-star 
radius ratio and corresponding transmission spectrum for 
near-grazing transits. 

It is difficult to diagnose the origin of the disagreement 
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between ATL AS and PHOENI X stellar atmosphere models 
for GJ 436; iKurucd (|2005[ ) note that the ATLAS mod- 
els should not be used for stellar effective temperatures 
below 3500 K, as they do not include important low- 
temperature opacity sources such as TiO and VO. How- 
ever, these molecules primarily affect the star's visible 
and near-infr ared spectra, and at 3 500 K they are still rel- 
atively weak ([Gushing et al.|[2005t ). Both disk- integrated 
and intensity spectra for the ATLAS models in this tem- 
perature range are featureless longward of 2.4 fim, with 
the exception of the GO band between 4.3 — 5.0 /im, 
whereas PHOENIX spectra show also clear molecular band 
structures, mainly due to H2O and OH, between 2.5 — 3.6 
/um and 6.5 — 8.0 /zm, with corresponding increases in the 
amount of limb darkening in these bands. The presence 
of the GO band in both model sets would appear to ex- 
plain the relatively good agreement in limb-darkening 
proles for the 4.5 /im Spitzer bandpass, but we were un- 
able to determine the reason behind the missing mid-IR 
H2O absorption features in the ATLAS models, which in- 
corporate the st rongest water lines CK urucz 1999) from 
the Ames list of [Partridge &: Schwenkel (|1997|). It is pos- 
sible that the spherical geometry used in the PHOENIX 
models (ATLAS models use a plane-parallel geometry) 
may also affect the resul t ing limb- darkening profile s 
(|Orosz k HauschildH [2000t IGlaret fc Hauschildtl [200l . 
but we find that PHOENIX models computed with a 
planet-parallel geometry show nearly identical limb dark- 
ening, with the exception of an exponential decline in 
the optically thin limb. We conclude that the differing 
opacities in the 3.6 and 8.0 /xm bands appear to be the 
most likely explanation for the disagreement between the 
limb-darkening profiles at these wavelengths. In this case 
the change in the value indicates the differences be- 
tween the two models are not statistically significant for 
this data set; near-IR grism spectr oscopy of transits o f 
GJ 436b, such as those obtained bv lPont et all (|2008bD . 
might help to better distinguish between these models. 

2.3.2. Error Analysis 

We calculate uncertainties for our best-fit transit pa- 
rameters using a M arkov Ghain Monte Garlo (MG MG) 
fit (see, for example iFoTdl 12001 IWhln et al.ll2007b( l with 
a total of 6 X 10^ steps, fourteen independent chains, and 
53 free parameters. Free parameters in the fits include; 
a/i?*, i, eight individual estimates of Rp/Ri,, eight tran- 
sit times, eight constants, a linear function of time and 
linear and quadratic terms in x and y for each of the 3.6 
and 4.5 ^m transits (20 variables total), the amplitude 
C2 and decay time C3 from Eq. [2] for the exponential fits 
to three 8.0 /im transits (6 variables total), and a linear 
function of time for the other 8.0 ^m visit. We assume 
a constant error for the points in each individual tran- 
sit light curve, defined as the the uncertainty needed to 
produce a reduced equal to one for the best-fit transit 
solution. 

We initialize each chain at a position determined by 
randomly perturbing the best-fit parameters from our 
Levenberg-Marquardt minimization. After running the 
chain, we search for the point where the value first 
falls below the median of all the values in the chain 
(i.e. where the code had first found the optimal fit), 
and discard all steps up to that point. We calculate the 
uncertainty in each parameter as the symmetric range 



Table 5 

Global System Parameters 



Parameter 



Value 



(d)b 
(d) 



TVansit Parameters 

in 

Duration T14 
T12 (Si T34)'^ 

b 

Rp (RisT 

T, (BJD) 

p id) 

Secondary Eclipse Parameters 
8 fim depth 
71 . d 

bright 
Duration T14 (d)° 
T12 (~ T34)<= (d) 
Orbital phase 
e cos{i.j) 
Tc(0) (BJD) 
P(d) 



86.699inro 

0.08311 ± 0.00026 
0.04227 ± 0.00016 
0.01044 ± 0.00014 
0.8521 ±0.0021 
0.0287 ± 0.0003 
0.437 ±0.005 
3.96 ± 0.05 
2454865.083208 ± 0.000042 
2.6438979 ± 0.0000003 



0.0452% ± 0.0027% 
740 ± 16 K 
0.04347 
0.00700 
0.58672 ± 0.00017 
0.13775 ±0.00027 
2454866.63444 ± 0.00082 
2.6438944 ± 0.0000071 



^ Calculated from the error-weighted average of the four 
8 /im planet-star radius ratio; this value was used for secondary 
eclipse fits. 

^ The transit duration T14 is defined as the time from first to 
fourth contact (i.e., the start of ingress to the end of egress). T12 
is the length of ingress, which is equal to the egress length in the 
limit of a circular orbit. Our best-fit transit ingress and egress 
lengths differ by less than 3 s. 

^ These parameters incorporate the stellar mass estimate of 
0.452 ± 0.013 Mq from lTorrel 120071) . 

^ Brightness temperature is defined as the temperature required 
to match the observed planet-star flux ratio in the 8 /im Spitzer 
band assuming that the planet radiates as a blackbody and using a 
Phoenix stellar atmosphere model (T^ff — 3585 K, log((7)— 4.843) 
for the star. 

^ The secondary eclipse duration and the length of ingress/egress 
were fixed in our fits. 

about the median containing 68% of the points in the 
distribution, except for the inclination and a/R^,, which 
we allow to have asymmetric error bars spanning 34% 
of the points above and below the median, respectively. 
The distribution of values was very close to symmetric 
for all other parameters, and there did not appear to be 
any strong correlations between variables. As a check 
we also carried out a residual permutat ion error analysis 
(|Gillon et al.l I2007bl : IWinn et al.l[2008h . which is sensi- 
tive to correlated noise in the light curve, on each in- 
dividual transit. At the start of each new permutation, 
we randomly drew values for the inclination and a/R^, 
from the simultaneous MGMG distribution and then fit 
for the corresponding best-fit values for the transit time 
and Rpl Ri, in that step. This ensures that our result- 
ing error distributions for individual transit times and 
Rp/ R-k values also take into account the uncertainties in 
the best-fit values for the inclination and a/Ri,. In each 
case where both a MGMG and residual permutation un- 
certainty are available for a given parameter we use the 
higher of the two values. We find that the MGMG fits 
generally produce larger uncertainties for the 8 fj-m ob- 
servations, whereas for 3.6 and 4.5 /xm data sets, which 
have higher levels of correlated noise, the residual per- 
mutation uncertainties are typically 50% larger than the 
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MCMC errors. 

2.3.3. Secondary Eclipse Fits 

We fit the secondary eclipses individually using the 
best-fit values for inclination and a/i?* from our tran- 
sit fits but allowing the eclipse depths and times to vary 
freely. Figure |4] shows the final binned data from these 
fits with the best fit normalizations for the detector ramp 
in each channel overplotted, and Figure [5] shows the 
binned data once these trends are removed, with best- 
fit eclipse curves overplotted. Best-fit parameters for 
individual eclipses are given in Table [3l and the error- 
weighted average (i.e., weights equal to l/cr^) of these 
eclipse depths is reported in Table [Sj We find that fixing 
the time of eclipse to a constant value, defined here as 
the best-fit orbital phase, does not significantly change 
our best-fit eclipse depths, nor does it reduce the uncer- 
tainties in our measurement of those depths. We calcu- 
late the uncertainties on individual eclipses using both a 
MCMC analysis and a residual permutation error analy- 
sis, again taking the higher of the two values as the final 
uncertainty for each parameter. 

3. RESULTS 

3.1. Orbital Ephemeris and Limits on Timing 
Variations 

We fit the transit times given in Ta ble [H to- 
gethcr with the transit times published in Pont et al.l 
j2008a); Bean & Scifahrt (2008); Coughhn ct al. (2008,)- 
lAlonso et a l. ( 2008 1 : IShporer et all (|2009( ) : ICaceres' et ahl 
(|2009D : iBaUard et air(r2010aD . with the following equa- 
tion, 

r,(n)=r,(0)+nxP (3) 

where Tc is the predicted transit time as a function of 
the number of transits elapsed since T'c(O) and P is the 
orbital period. We find that Tc = 2454865.083208 ± 
0.000042 BJD and P = 2.6438979 ± 0.0000003 days. As 
demonstrated by Fig, O the 34 published transit times 
appear to be markedly inconsistent with a constant or- 
bital period, with the most statistically significant out- 
liers (6.2 and 7.1a, respectively), occurring during the 
sequence of eight transits observe d by the EPOXI mis - 
sion between UT 2008 May 5-29 (|Ballard et al.ll2010af ). 
The most significant deviations in the Spitzer transit 
data presented here occur during the last three visits 
(UT 2009 Jan 28 - Feb 2), and range between -3.1 and 
-|-3.5(T in significance. Given the size of these discrep- 
ancies, it is perhaps not surprising that the reduced 
value for the linear fit to Eq. [3] is 6.8 (total of 216.4, 
34 points, two degrees of freedom). It is unlikely that the 
observed deviations could be explained by perturbations 
from a previously unknown second planet in the system, 
as the measured transit times shift by as much as several 
minutes on time scales of only a few days (i.e., a single 
planet orbit). As we discuss in more detail in §4.1.3) we 
believe the presence of occulted star spots in a subset of 
the transit light curves is the most likely explanation for 
the observed deviations. 

We carry out a similar fit to the secondary eclipse 
times given in Table |3l along with the additional sec- 
ondary eclipse times reported in .Stevenson ct al.. (2010i ). 
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Figure 8. Observed minus calculated transit times using the new 
best-fit ephemeris. The dashed lines indicate the ztltr uncertainty 
in the predicted transit times, with a solid line at O — C equal to 
zero. Spitzer measurements from this paper are plotted as filled 
circles, and previously published observations are shown as filled 
stars. The color of the points denotes the wavelength of the obser- 
vations (blue for visible, red for IR). Moving f rom left to right, tran- 
sits between 200-300 BJD-2454000 are from 'Shporer et al. ( 2003) 
and Caceres et al. (2009), transits between 400-500 BJD-2454000 
with small uncertainties are from 'Pont et al. ( 2008a) and those 
with large uncertainties are from Bean & Seifahrt (2008 ). Be- 
tween 530-620 BJD-2454000, observations arc from Co ughlin et al.l 
(|2008), Alonso et al. (2008), and Ballard et al. (2010^). Visible- 
light transit observations typically show larger timing variations 
than the IR observations, indicating that spot occultations may be 
responsible for the apparent timing variations. 

and find that Tc{Q) ^ 2454866.63444± 0.00082 BJD and 
P = 2.6438944 ± 0.0000071 days. This period is consis- 
tent with the best-fit transit period to better than la, 
and we therefore conclude that there is no evidence for 
orbital precession in this system. We also see no evi- 
dence for statistically significant variations in the sec- 
ondary eclipse times (see Fig. [3]), as would be expected 
if the shifted transit times were due to occulted spots, 
but our measurements are not precise enough to rule out 
timing variations of the same magnitude as those ob- 
served in the transit data. If we fix the orbital period 
to the value from the transit fits and subtract the 28 s 
light travel time delay for this system (Locb 2005'), we 
find that the secondary eclipses occur at an orbital phase 
of 0. 58672 ± 0.00017. consis tent with the best-fit phase 
from iStevenson et al.l (|2010D . 

We can use the offset in the best-fit secondary eclipse 
time to calculate a new estimate for ecos(a;). We find 
that the secondary eclipse occurs 330.18 ± 0.67 minutes 
later on average than the predicted time for a circular or- 
bit, including the correction for the light travel time. We 
can convert this to e co s (a;) u sing the expression reported 
in Eq. 19 of iPal et all (|2010D . Note that this expression 
is more accurate than the commonly used approxima- 
tion of ecos(w) ~ ^ fe.g.. iCharbonneau et al.l [20051 : 
iDeming et al.| [2b05). where 5t is the delay in the mea- 
sured secondary eclipse time and P is the planet's orbital 
period. We find that using the less accurate approxima- 
tion gives ecos(a;) = 0.13622 ± 0.00026, while the equa- 
tion from Pal et al. yields e cos(w) ^ 0.13754±0 . 00027, a 
4g di fference in this case (see also, ISternel[T940t Ide KortI 
Il954f ). If we take the best-fit longitude of pericenter from 
the radial velocity fits, 334° ± 10°, we find an orbital ec- 
centricity equal to 0.153 ±0.014. This is consistent with 
the current best-fit orbital eccentricity from radial veloc- 
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Figure 9. Observed minus calculated secondary eclipse times us- 
ing the best-fit period from the transit fits and allowing the phase 
of the secondary eclipse to vary freely. Filled circles are eclipse 
times reported in this paper, and filled stars are additional 3.6 
and 4.5 /im eclipse times from [Stevenson et al.l (12010 ^ The solid 
line indicates the best-fit phase, with itlcr uncertainties plotted as 
dashed lines. 



Table 6 

a/R* and Inclination Values from Independent Transit 
Fits 



UT Date A (^tm) 




i 


n 


a/R 




UT 2007 Jun 29 


8.0 


86.68 


± 


0.12 


14.11 




0.35 


UT 2008 Jul 14 


8.0 


86.54 


± 


0.12 


13.67 




0.36 


UT 2009 Jan 9 


3.6 


86.76 


± 


0.07 


14.40 


± 


0.22 


UT 2009 Jan 17 


4.5 


86.85 


± 


0.10 


14.60 




0.32 


UT 2009 Jan 25 


8.0 


86.70 


± 


0.14 


14.19 


± 


0.43 


UT 2009 Jan 28 


3.6 


86.67 


± 


0.07 


13.96 




0.20 


UT 2009 Jan 30 


4.5 


86.58 


± 


0.10 


13.76 




0.28 


UT 2009 Feb 2 


8.0 


86.80 


± 


0.14 


14.49 


± 


0.44 



ity data alone, e = 0.145 ± 0.017 (A. Howard, personal 
communication, 2010). 

3.2. System Parameters from Transit Fits 

In this work we examine two transits obtained at 
3.6 //m, two transits at 4.5 /im, and four transits at 
8.0 fjim. We carry out two sets of transit fits, one where 
the ratio of the orbital semi-major axis to the stellar 
radius a/i?* and the orbital inclination i are allowed 
to vary freely, and the other where they have a sin- 
gle common value for all visits. In all cases we allow 
the planet-star radius ratio and best-fit tran- 

sit times to vary independently for each visit. In fits 
where a/Ri, and i are allowed to vary individually we 
find no evidence for statistically significant variations in 
either of these parameters (see Table |6]) , and we there- 
fore proceed assuming that these parameters have a sin- 
gle common value in our subsequent analysis. Our best- 
fit values for i, a /Rj,, and Rp/R^, are consistent with 
those reported bv I Ballard et al.l ()2010aD to better than 
la, and the impact parameter b and transit duration 
T = Ti4 - Ti2 = 0.0318 ± 0.0007 days that we derive 
from our fits are si milarly consistent with the value re- 
ported bv lPont et a l. (2008b). 

Although the best-fit orbital inclination and a/Ri, ap- 
pear to be consistent with a constant value over the 
approximately two year period spanned by our obser- 
vations, we do see evidence for statistically significant 
differences in the transit depths within the same Spitzer 



bandpass (see Fig. [TU)) . We would expect to see the 
transit depth vary with wavelength due to absorption 
from the planet's atmosphere, but this signal should re- 
main constant from epoch to epoch for observations in 
the same bandpass. If we compare individual visits in 
a given bandpass, we find that the two 3.6 ^m radius 
ratios, measured on UT 2009 Jan 9 and 28, are incon- 
sistent at the 4.7(7 level. The two 4.5 /im radius ratios, 
measured on UT 2009 Jan 17 and 30, differ by 2.9ct. The 
four 8 /im transits, measured on UT 2007 Jun 29, UT 
2008 Jul 14, UT 2009 Jan 28, and UT 2009 Feb 2, differ 
from the error-weighted average by 0.2cr, I.Oct, 1.5ct, and 
2. Oct, respectively. These offsets are still present in the 
fits where the inclination and a/Ri, are allowed to vary 
individually, indicating that the discrepancy cannot be 
due to a change in these two parameters. 



4.1. 



4. DISCUSSION 

Transit Depth Variations 

In the sections below we consider three possible expla- 
nations for the observed depth variations: first, that the 
effective radius of the planet is varying in time, second, 
that residual correlated noise in the data affected the 
best-fit transit solutions, and third, that spots or other 
stellar activity produced apparent depth variations. 



4.1.1. A Time- Varying Radius for the Planet 

We first consider the possibility that the radius of the 
planet is changing in time, either due to thermal expan- 
sion of the atmosphere or the presence of a variable cloud 
layer at sub-mbar pressures. We require a change in ra- 
dius of approximately 4% in order to match both of the 
measured 3.6 /im transit depths; if this change is due 
to thermal expansion, we can estimate the energy input 
required using simple scale arguments. 

The effective change in the planet's radius due to heat- 
ing of the atmosphere depends on both the amount of 
heating and the range in pressures over which this heat- 
ing t akes place. We use the secondary eclipse depths in 
i i4.3l to place an upper limit on the allowed change in 
temperature at the level of the mid-IR photosphere, and 
then calculate the corresponding range in pressure that 
must be heated by this amount in order to increase the 
radius of the planet by 4%. If we assume a hydrogen at- 
mosphere with a baseline temperature of 700 K, we find 
a corresponding scale height of approximately 240 km, 
where the scale height is defined as H = —, T is the 

t> ^lg ^ 

temperature of the atmosphere, ^ is the mean molecu- 
lar weight, and g is the surface gravity. We know from 
the secondary eclipse observations described in i )4.3l that 
the temperature of the planet's dayside atmosphere must 
change by less than 30%, which would correspond to an 
upper limit of 100 km on corresponding changes in the 
planet's scale height. 

In order to calculate the required energy input to pro- 
duce the observed change in radius, we must first deter- 
mine the range of pressures affected by this heating. We 
model the planet as an interior region with a constant 
temperature, surrounded by an outer envelope that ex- 
pands and contracts freely with changing temperature. 
We set the upper boundary on this region equal to 50 
mbar, corresponding to the approximate location of the 
T = 1 surface in the mid-infrared. As illustrated in Fig. 
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[Til opaque clouds at this pressure suppress but do not 
entirely remove absorption features in the planet's trans- 
mission spectrum at these wavelengths, making this a 
reasonable estimate for the location of the r = 1 sur- 
face. We assume that when the planet is heated the 
scale height changes by 100 km, which requires the lower 
boundary of the heated region to be located at a pressure 
of approximately 1 bar in order to produce a 1% expan- 
sion in radius. If we then calculate the change in the 
planet's gravitational energy corresponding to this ex- 
pansion, we find that an energy input of approximately 
10^^ J is required. Repeating this calculation for a 4% 
increase in radius, we find a lower boundary at 8,000 
bars and a corresponding energy input of 10^^ J. The 
insolation received by the planet is 10^" W, which gives 
an energy budget of 10^^ J per orbit. When we examine 
Fig. [10] we find that that the observed change in radius 
occurs primarily between the third and fourth visits (UT 
2009 Jan 25-28). This would require an energy input as 
much as 10^ times higher than the total insolation over 
this epoch, which is clearly unphysical. 

One alternative explanation for the observed change in 
radius would be to invoke the presence of intermittent, 
high-altitude clouds. Such clouds could produce a change 
in the apparent radius of the planet across multiple bands 
without requiring any actual heating or cooling of the 
atmosphere. In this picture, smaller radii for the planet 
would correspond to the cloud-free state, while larger 
radii would require the presence of an additional cloud 
layer. A change of 4% in apparent radius would require 
the clouds to form at a pressure approximately 100 times 
lower than the location of the nominal cloud-free radius. 
In §4.21 we find that the average pressure of the r = 1 
surface for the nominal methane-poor (green) model be- 
tween 3 — 10 /im is 40 mbar, indicating that the clouds 
would have to extend to 0.4 mbar to explain the largest 
measured 3.6 /im radius for the planet. This conclusion 
is reasonably independent of our assumed composition, 
as the average r = 1 surface for the methane-rich (blue) 
model is located at 30 mbar. Gravitational settling would 
presumably pose a challenge for cloud layers at sub-mbar 
levels, but vigorous updrafting of condensate particles 
might compensate for this effect. The broadband nature 
of the data presented here make it difficult to directly test 
this hypothesis; we therefore recommend the acquisition 
of high signal-to-noise, near-infrared grism spectroscopy 
over multiple transits in order to resolve this issue. A 0.5 
mbar cloud layer would lead to a near-featureless trans- 
mission spectrum whereas a lower cloud layer would still 
exhibit many of the same absorption features as a cloud- 
free atmosphere. Such a data set would also allow us 
to test the theory, outlined in > j4.1.31 that the observed 
transit depth variations are due to the occultation of 
regions of non-uniform brightness on the surface of the 
star, as these regions should also produce a wavelength- 
dependent effect. 

4.1.2. Poorly Corrected Systematics 

It is possible that poorly corrected instrument effects, 
such as the intrapixel sensitivity variations at 3.6 and 
4.5 micron, or the detector ramp at 8.0 /xm, might lead 
to variations in the measured transit depth. Because 
there is complete overlap between the positions spanned 
by the star in the in-eclipse and out-of-eclipse data for all 
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Figure 10. Best-fit transit (upper panel) and secondary eclipse 
(lower panel) depths as a function of time. Average transit /eclipse 
depths are shown as a dashed line in each plot. 3.6 ixra observations 
are denoted with stars, 4.5 /im observations are denoted with tri- 
angles, and 8.0 /xm observations are marked with solid circles. The 
most recent three transits are systematically high when compared 
to earlier transit visits; this appears to coinci de w ith a decrease in 
the total visible light flux from the star (Fig. Illll . suggesting that 
the fractional spot coverage on the star's visible face was increasing 
in time during the later part of our observations. 



3.6 and 4.5 /.tm visits, fits that inadequately describe the 
pixel response as a function of position should fail equally 
for both sections of the light curve. The UT 2009 Jan 
30 transit serves as an example of imperfectly removed 
detector effects, as the residuals display a sawtooth sig- 
nal with a shape and timescale similar to the original 
intrapixel sensitivity variations (see ^2.1\ for a more de- 
tailed discussion of this light curve). Conversely, it is 
much more difficult to explain the 3.6 micron transit on 
UT 2009 Jan. 28 with this scenario, as there appears to 
be a large dip in residuals during ingress, but when the 
star spans the same pixels in the out-of-eclipse data we 
see no comparable deviations. 

In this section we consider an alternate decorrela- 
tion function that better accounts for small-scale varia- 
tions in the in t rapixel sensitivity function as discussed in 
[Ballard et al.l (|2010b[ ). Following the discussion in Bal- 
lard et al., we describe the intrapixel sensitivity varia- 
tions using a position- weighted average of the time series 
after the best-fit transit function and linear function of 
time from the position fits described above has been di- 
vided out. Unlike Eq. [l] this formalism does not assume 
a functional form for the intrapixel sensitivity variations 
and therefore should in principle produce an unbiased 
correction for these variations. We calculate the weight- 
ing function as follows: 
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Eexp(-^^ 



2^ 



(4) 



where Xi and yi are the x and y positions of the i^^ frame, 
Xj and j/j are the x and y positions for the rest of the time 
series. We optimize our choice of Ux and Uy to produce 
the smallest possible scatter in the final time series when 
we fix the transit light curve to the best-fit solutions 
listed in Table [2] We find that the preferred values range 
between 0.0053 - 0.0120 pixels in and 0.0024 - 0.0045 
pixels in Uy for the four 3.6 and 4.5 /im transits examined 
here. For ease of computation we bin our time series in 
intervals corresponding to one point per original set of 64 
images (in some instances there are less than 64 images 
in a given bin after removing outliers) and iteratively 
calculate the weighting function and the linear function 
of time plus transit fits until we converge to a consistent 
solution. 

Once we have a final solution we calculate the weight- 
ing function for the unbinned data and carry out a final 
fit for the transit function to determine our best-fit tran- 
sit depth. In this case we fix the inclination and a/ Ri, 
to their best-fit values from the simultaneous fits to all 
transits described in ^2.'i\ which allows us to fit each 
transit individually using the weighting function while 
still preserving the constraints imposed in a simultane- 
ous fit. We find that in all cases we obtain transit depths 
and times that are consistent with the values from our 
fits using Equation [1] with a standard deviation that is 
comparable or slightly worse than that achieved with our 
polynomial fits. 

We also carried out a second set of fits in which we 
derived our corrections for the intrapixel sensitivity vari- 
ations using only the out-of-transit data, and found that 
our best-fit planet-star radius ratios changed by less than 
0.4cr in all cases. Because the star samples the same re- 
gions of the pixel in both the in-transit and out-of-transit 
data, it is possible to obtain an equivalently good correc- 
tion for the intrapixel sensitivity variations using only the 
out-of-transit points. Conversely, this means that poor 
corrections for this effect should produce equally large de- 
viations in both the in-transit and out-of-transit regions 
of the light curve. As we will discuss in the following sec- 
tion, we find that the residuals for the deepest transits in 
these two bands have a significantly higher RMS in tran- 
sit than out of transit. This behavior is inconsistent with 
our expectations for poorly corrected instrument effects, 
and we therefore conclude that it is unlikely that these 
effects are responsible for the discrepant transit depths 
measured at 3.6 and 4.5 /zm. 

At 8.0 /im we fit the data with a single or dou- 
ble exponential function to describe the s moothly vary- 
ing detector ramp. In lAgol et all (|20T?! ) we conclude 
that this functional form avoids correlations between 
the slope of the ramp and the measured transit or 
eclipse depth; however we check this assertion using 
our 8 /Ltm data as well. For our 8 /im transit fits 
we find that the exponential term has a coefficient of 
[0.00156, 0.00000, 0.00288, 0.00299], corresponding to 
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Figure 11. The upper panel plots APT measurements of vari- 
ation in the averaged Stromgren b+y band fluxes (filled circles) 
obtained for GJ 436 between UT 2008 Nov 30 and UT 2009 May 
29, where the epoch of our 2009 Spitzer observations is denoted 
by the grey shaded region. These bandpasses are sensitive to the 
rotation-modulated flux of the star, which we flnd has a best-flt 
period of 57 days during this epoch. We overplot a red curve show- 
ing our best sine-curve fit for the spot modulation, together with 
a quadratic function of time to describe the evolution of the spot 
coverage on longer time scales. The lower panel shows the mea- 
sured values for G J 436 from th e Keck HIRES instrument 
during this same period lllsaacson fc Fi scher 2010), with a best-fit 
sine -I- quadratic function overplotted in red. Error bars for both 
panels are set equal to the standard deviation of the residuals. Al- 
though the best-fit period of 57 days for the Sjjpj data during this 
epoch is only marginally significant, the S^j^ values appear to be 
anti-correlated with the flux variations. This is consistent with a 
model in which increased magnetic activity is associated with the 
presence of spots or other dark regions on the surface of the star. 

planet-Star radius ratios of [0.08234, 0.08162, 0.08138, 
0.08336], where we have set the amplitude of the ex- 
ponential term to zero for the transit occurring in the 
middle of our 70-hour phase curve observation. For the 
eleven secondary eclipse observations we find coefficients 
of [0.00645, 0.00627, 0.00194, 0.00433, 0.00534, 0.00140, 
0.00000, 0.00359, 0.00321, 0.00353, 0.00262], correspond- 
ing to eclipse depths of [0.0552, 0.0507, 0.0395, 0.0495, 
0.0367, 0.0523, 0.0421, 0.0386, 0.0491, 0.0397, 0.0441], 
respectively, where we have set the exponential coeffi- 
cient to zero for the secondary eclipse at the end of the 
phase curve observation. We find no evidence for any 
correlation between the slope of the exponential function 
and the measured transit or secondary eclipse depths. As 
an additional check we also confirm that there is no cor- 
relation between these depths and cither the measured 
sky background or the total stellar flux given in Table [TJ 

4.1.3. Stellar Variability 

The presence of spots or faculae on the visible face 
of the star can have two distinct effects on the mea- 
sured light curve for a transiting planet. Non-occulted 
spots on the visible face of the star reduce the star's 
total flux, increasing the measured transit depth, while 
spots occulted by the planet cause a small positive devi- 
ation in the light curve with a time scale proportional to 
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Figure 12. Averaged Stromgren b band fluxes vs. ^hK 
GJ 436 over a period of seven years. Because tiie sampling for 
the 5jjj^ measurements is much lower than for the flux measure- 
ments in most seasons, for the purposes of this plot the b fluxes 
are defined as the average of all b band measurements taken within 
one day of the Sjjj^ measurement. We also average 5jjj^ values 
take n on the same night, scaling the flux and ^jj]^ error bars from 
Fig. Illl bv the square root of the number of measurements in each 
bin. We find that using b photometry instead of the averaged b+y 
fluxes results in increased scatter but also strengthens the observed 
correlation. 

the p hysical size of the occulted spot (e.g., IRabus et al.l 
l2009f l. and occulted faculae would have the opposite ef- 
fect. The early K dwarf HD 189733 (Toff = 5100 K) 
is perhaps the best-studied example of an active star 
with a tra nsiting hot Jupiter (e.g.. iBakos et al.l 120061 : 
iPont et al.l l2008at iDesert et all |2011aD, but the late G 
dwarf CoRoT-2 (Toff = 5600 K) also exhibits a high level 
of spot activity that may have resulte d in early over- 
estimates of its planet's inflated radius (jGuillot fc Havell 
120111) . This problem is likely to be even more common for 
M dwarfs, and in fact several instances of occulted spots 
were reported in transit light curves for the super -Earth 
GJ 1214b, which orbits a 3000 K primary (Cart er et al.l 
[Ml IBerta et al.l[201ll : IKundurthv eIaLl[20ilf l. 

Although it is important to correct for these effects in 
any transit fit, it is particularly crucial when comparing 
non-simultaneous, multi-wavelength transit observations 
such as the ones described in this paper, which require 
a relative precision of better than a part in 10~^ in the 
measured transit depths. We evaluate the likely impact 
of GJ 436's activity on the measured transit depths us- 
ing several complementary approaches. First, we esti- 
mate the average activity level on GJ 436 by measuring 
the amount of emission in the cores of the Ca II H & 
K lines; in Knutson et al. (2010) we determined that GJ 
436 had an average 5hk of 0.620. Ilsaacson fc Fisclieil 
(pOlQ) found that other stars in the California Planet 
Search database with similar B — V colors have ^hk val- 
ues ranging between 0.5 — 2.0, indicating that GJ 436b 
is rela tively quiet for its spectral type. iDemorv et al.l 
(|2007| ) report that this star's rotation period is greater 
than 40 days, consiste nt with upper limits on v sin i of 
1 km/s (| Jenkins et al.l 12009'). also suggesting that it is 
likely to be relatively old and correspondingly quiet. 
The upper limit of 3 km/s on usini from spectroscopy 
(|Butler et al.ll200"l) is also consistent with an inclined or 
pole-on viewing geometry, although it is not required as 



long as the star's rotation period is longer than 7 days. 

Rather than relying on these indirect measures of ac- 
tivity, we can also directly measure the amplitude of the 
star's rotation-modulated flux variations using visible- 
light ground-based observations. We obtained observa- 
tions of GJ 436 in Stromgren b and y filters over a span of 
approximately six months surrounding our 2009 Spitzer 
transit and secondary eclipse observations from an on- 
going monitoring program carried out with the T12 0.8 
m AP T at Fairborn Observato ry in sou thern Arizona 
(|Henrv 1999; Eaton et al. 2003; Henrv fc W inn 2008). In 
these observations the telescope nodded between GJ 436 
and three comparison stars of comparable or greater 
brightness, which were then used to correct for the ef- 
fects of variable seeing and airmass. We find that during 
the period between UT 2009 Jan 9 - Feb 4, when a ma- 
jority of our transit data was obtained, the star varied 
in flux by less than a few mmag in visible light (Figure 
[TT|) . We carry out a similar check for variability in the 
infrared using the fifteen 8 /zm flux estimates listed in Ta- 
ble [1] which we find have a standard deviation of 0.07%. 
Both of these measurements indicate that the star is very 
nearly constant in flux in both visible and infrared light, 
and we can therefore rule out non-occulted spots as the 
cause of the observed transit depth variations. 

We also use these same data to search for periodici- 
ties corresponding to GJ 436b's rotation period. If we 
fit the combined b and y band fluxes with a sine func- 
tion plus a quadratic function of time as shown in Fig. 
[TTl we find a best-fit period of 56.5 days. We calculate 
a Lomb-Scargle periodogram (lLomblll976l: IScarglel[T98l) 
for these data and find that this period has a false alarm 
probability of only 2%, which we determine using a boot- 
strap Monte Carlo analysis. We find a nearly identical 
best-fit period of 56.6 days in the 5 hk values measured 
with Keck HIRES during this epoch ([Isaacson fc Fischen 
l2010f) . but with a much higher false alarm probability 
of approximately 20%. We also examine the correlation 
between the measured b fluxes and S'hk values over the 
six- year period in which both were available (Fig. I12[) . 
and find that these parameters are negatively correlated. 
Taken together, these data indicate that the small ob- 
served variations in GJ 436's visible-light fiuxes are likely 
connected with the presence of regions of increased mag- 
netic activity on the visible face of the star. 

Although such low-amplitude fiux variations generally 
indicate that a star has relatively few spots, there are 
two important exceptions. First, if the spots are uni- 
formly distributed in longitude, it is theoretically possi- 
ble to have a star with significant spot coverage and an 
effectively constant fiux. It would not be surprising if the 
occurrence rate and distribution of spots was different for 
M stars than for G or K stars, but in GJ 436b's case the 
lack of any fiux variations larger than a few mmag would 
seem to place a strong limit on the allowed spot distribu- 
tions. We can quantify this limit if we assume that the 
deviation of approximately 0.08% in the first part of the 
3.6 /im transit light curve from UT 2009 Jan 28 shown 
in Fig. [2] is due to the occultation of a bright region on 
the star. This region must have a surface intensity that is 
12% brighter than the rest of the star in order to produce 
the observed deviation. If we compare PHOENIX models 
with varying effective temperatures integrated over this 
band, we find that the star's temperature must increase 
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by approximately 200 K in the affected region in order 
to match this surface intensity. We know that the total 
rotational modulation in the star's visible-light flux must 
remain below 0.1%, and we estimate that an increase of 
12% in the 3.6 fim surface intensity should produce an 
increase of approximately 65% in the Stromgren {b+y)/2 
band. In this case the fractional area covered by active 
regions on the star must vary by less than 0.15% from 
the most active to the least active hemisphere. Of course, 
it is possible that the stellar atmosphere models do not 
provide an accurate match for the spectra of these active 
regions; if we instead use the measured 3.6 /im flux con- 
trast of 12%, we find a more conservative limit of 1% on 
variations in the area affected by stellar activity. 

A second, more plausible scenario involves tilting the 
rotation axis of star so that we are viewing it closer to 
pole-on, which would effectively suppress the amplitude 
of rotational flux variations regardless of spot coverage. 
If we assume that the star's spin axis is randomly ori- 
ented with respect to our line of sight, the probability 
that it will fall within 45° of a pole-on view is 30%. In this 
scenario the star could be highly spotted, allowing for fre- 
quent occultations of spots by the planet, while still dis- 
playing a small rotational flux modulation. This scenario 
would require the planet's orbit to be misaligned with re- 
spect to the star's rotation axis, but such misalignments 
are commo nly seen in other transiting planet systems 
(jWinn et al . 2010a). Although the Rossiter-McLaughlin 
effect has never bee n successfully measured for GJ 436b, 
I Winn et al.l (|2010bD find that the Neptune-mass planet 
HAT-P-llb, which is perhaps the best analogue to the 

4-26° 

GJ 436 system, has a sky-projected obliquity of 103° o 

indicating that this system is significantly misaligned. If 
most close-in planets start out misaligned and are then 
gradually brought into alignment through tidal interac- 
tions w ith their host star as proposed by Winn et al. 
(|2010af) . the fact that HAT-P-llb stiU maintains both 
a non-zero orbital eccentricity and a significant misalign- 
ment would seem to suggest that the same could also be 
true for GJ 436b. 

If we proceed with the hypothesis that GJ 436 is both 
spotty and tilted with respect to our line of sight, we can 
then search for evidence of occulted spots in the light 
curves with discrepant transit depths. We first compare 
the relative standard deviations of the in-transit (crin) 
and out-of-transit (cin) residuals plotted in Fig. [2j 



Cout 

We list the measured values of drei for all eight transit 
observations in Table [T] Both the 3.6 /im transit on UT 
2009 Jan 28 and the 4.5 fim transit on UT 2009 Jan 30 
appear to have inflated values of arei, as would be ex- 
pected if the planet occulted active regions on the star 
during these visits. We can quantify the statistical signif- 
icance of the measured cTj-ei values if we assume that both 
the in-transit and out-of-transit points are drawn from 
the same underlying Gaussian distribution, and then ask 
how many times in a sample of 100,000 random trials 
we measure a value of cr-rei greater than or equal to the 
value calculated directly from our observations. In each 
trial we generate two synthetic data sets, each with the 



Table 7 

A Comparison of the In- Transit vs. Out-of- Transit Standard 
Deviations 



UT Date 


A (/im) 




Aout'' 


CTrel 


P(o-rel) 


UTlbZTlTlGd DtttOj 














UT 2007 Jun 29 


8.0 


7,924 


17,956 


-1 


.3% 


0.92 


UT 2008 Jul 14 


8.0 


7,895 


25,012 


-1-1 


.4% 


0.059 


UT 2009 Jan 9 


3.6 


25,482 


81,848 


-1-0, 


.2% 


0.36 


UT 2009 Jan 17 


4.5 


25,954 


82,212 


-0, 


.3% 


0.70 


UT 2009 Jan 25 


8.0 


7,910 


25,334 


+0. 


.1% 


0.44 


UT 2009 Jan 28 


3.6 


25,536 


82,238 


+1. 


.4% 


0.0023 


UT 2009 Jan 30 


4.5 


25,955 


62,334 


+1. 


.1% 


0.018 


UT 2009 Feb 2 


8.0 


7,890 


25,318 


+0. 


.1% 


0.45 


Binned Data 














UT 2007 Jun 29 


8.0 


126 


287 


-13, 


.6% 


0.97 


UT 2008 Jul 14 


8.0 


126 


399 


-4, 


.9% 


0.75 


UT 2009 Jan 9 


3.6 


411 


1311 


+3. 


.3% 


0.21 


UT 2009 Jan 17 


4.5 


412 


1310 


-3, 


.4% 


0.80 


UT 2009 Jan 25 


8.0 


126 


403 


+9. 


.7% 


0.093 


UT 2009 Jan 28 


3.6 


411 


1311 


-1-37, 


.5% 


1 X 10-6 


UT 2009 Jan 30 


4.5 


412 


989 


-1, 


.4% 


0.63 


UT 2009 Feb 2 


8.0 


126 


403 


+2. 


.9% 


0.34 



^ Number of in-transit and out-of-transit points. 

^ Probability ttiat the standard deviation of the in-transit data would 
be greater than the standard deviation of the out of transit data by an 
amount CTrei if both data sets are drawn from the same underlying Gaus- 
sian distribution. 

appropriate length corresponding to either the in-transit 
or out-of-transit measurements, and then calculate the 
standard deviation of each distribution and the corre- 
sponding value of (Trei- In the 3.6 micron transit obser- 
vation on UT 2009 Jan 9 there are 81,848 out-of-transit 
flux measurements and 25,482 in-transit flux measure- 
ments, and we find that over 100,000 trials, we obtain a 
value of cTrei greater than or equal to the measured value 
of 0.2% approximately 36% of the time. Repeating the 
same calculation for the 3.6 micron transit observed on 
Ut 2009 Jan 28, which has 82,238 out-of-transit points 
and 25,530 in-transit points, we obtain a-rei greater than 
or equal to the measured value of 1.4% only 0.23% of the 
time. We list the corresponding probabilities for all eight 
transits in Table [71 

We also repeat this same test with data that has been 
binned in sets of 64 images, corresponding to 10 s bins at 
3.6 and 4.5 um and 30 s bins at 8 um. This allows us to 
evaluate the relative contribution that correlated noise 
makes to the in-transit and out-of-transit variances, as 
the photon noise should be be reduced by a factor of 8 
in these bins (also see Fig. 3). In this case we carry 
out 1,000,000 random trials for each visit, as each simu- 
lated data set is much smaller and the computations are 
correspondingly fast. We find that for the binned Jan 9 
light curve there are 1311 points out of eclipse and 411 
points in eclipse. In this case CTrei is 3.3%, and we ob- 
tain values greater than or equal to this number in 21% 
of our random trials. Repeating this calculation for the 
UT 2009 Jan 28 visit, we find that the measured value of 
frei is 37% (i.e., a standard deviation that is 37% higher 
in eclipse than it is out of eclipse), with 1311 points out 
of eclipse and 411 points in eclipse. In our simulations 
assuming a single Gaussian probability distribution for 
both segments, this level of disagreement occurred only 
once in 10^ trials. We find that in all other visits, includ- 
ing the 4.5 micron transit observed on UT 2009 Jan 30, 
the binned data in and out of eclipse are consistent with 
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a single distribution. 

One consequence of a misalignment between the star's 
rotation axis and the planet's orbit is that the planet will 
not necessarily occult the same spot on successive tran- 
sits, as would be expected for a well-aligned system; we 
therefore consider each transit individually. Our analysis 
above indicates that the 3.6 fim transit on UT 2009 Jan 
28 displays a statistically significant increase in the stan- 
dard deviation of the in-transit data that is dominated by 
contributions from correlated noise on time scales greater 
than 30 s, as would be expected if the planet occulted 
an active region on the surface of the star. Although the 

4.5 fj,m transit from UT 2009 Jan 30 does not appear to 
display a similar increase, our imperfect correction for 
the intrapixel sensitivity variations in this visit means 
that we are less sensitive to variations in CTrei- We argue 
that even if the star's rotation axis and the planet's or- 
bit are misaligned, it is still likely that the planet would 
occult the same active region during both the UT 2009 
Jan 28 and Jan 30 visits, as the interval between these 
visits is much shorter than the star's approximately 50 
day rotation period. As we discuss later in this section, 
the fact that both visits display increased transit depths 
and shifted transit times provides additional support for 
this hypothesis. 

We also consider the possibility that the increased scat- 
ter in the in-transit residuals might be due to a change in 
the transit parameters, including the planet's radius, or- 
bital inclination, transit time, or a/ R^,, from one visit to 
the next. We test this hypothesis by taking the difference 
of the first and second visits in each bandpass from 2009 
and comparing the shape of the residual light curve to 
the differences we would expect due to changes in these 
parameters, which should be distinct from the deviations 
created by occulted star spots (Fig. [T3|) . Because we are 
directly differencing the two light curves, our results are 
independent of any assumptions about the shape of the 
transit light curve or the stellar limb-darkening. We in- 
spect the deviations in the residuals plotted in Fig. [13] 
and conclude that they do not appear to be well-matched 
by changes in the best-fit transit parameters, leaving oc- 
cultations of active regions on surface of the star as the 
most likely hypothesis. 

If the planet occults a spot it can also cause a shift in 
the best-fit transit times, particularly when the spot is 
near the edge of the star and is occulted during ingress 
or egress. Indeed, we see that the UT 2009 Jan 28 

3.6 fim appears to occur 31.4 ± 9.5 s early, while the 
4.5 iim Jan 30 visit occurs 34.4 ± 9.4 s late (see Fig. [5]) 
in the fits where we fix a/i?* and i to a single common 
value. As a test we repeated our fit to the 3.6 /im tran- 
sit excluding the first 1/3 of the transit light curve, and 
found that the best-fit transit time shifted forward by 
approximately 30 s. We would also expect that tran- 
sits observed in visible light, where the contrast between 
the spots and the star is more pronounced, would show 
proportionally larger timing deviations when the planet 
crosses a spot. As noted in tj3.11 the scatter in the mea- 
sured visible-light transit times is inconsistent with a 
constant period, and the amplitude of the visible-light 
deviations is on average larger than the deviations in 
the infrared. We should also see this same wavelength- 
dependence in the measured transit depths in Fig. llOi 
and indeed we find that the 3.6 /im transit depth changes 



by 7.8%, the 4.5 fim transit depth changes by 5.3%, and 
the 8.0 /zm transit depth changes by 4.9% during the pe- 
riod between UT 2009 Jan 9 and Feb 2. Lastly, we can 
examine the visible-light fiux measurements for GJ 436 
in Fig. [TT] and see that these two transits were obtained 
near a minimum in the star's flux, consistent with a rela- 
tive increase in the fractional spot coverage as compared 
to earlier epochs. The measured values for S'hk, a com- 
mon activity indicator, appear to be anti-correlated with 
the observed flux variations and reach a local maximum 
near this point. 

4.2. Atmospheric Transmission Spectrum 

In principle, the broadband transmission photometry 
of GJ 436b allows us to constrain the chemical composi- 
tion and temperature str ucture near the limb of the plan - 
etary atmosphere (e.g.. 'Mad husudhan fc Seageri 120091) . 
However, time variability in either the properties of the 
star or of the planet poses a significant challenge to an 
analysis in which we are comparing transit observations 
at different wavelengths obtained days or weeks apart. 
As discussed in i )4.1 11 we consider it unlikely that the 
discrepancies in the measured transit depths arc due to 
changes in the properties of the planet, but instead con- 
clude in H4.1.3l that the occultations of regions of nonuni- 
form brightness in a subset of the transits appear to be 
responsible for the observed depth variations. 

If we set aside those transits which we believe to be 
most strongly affected by stellar activity, including the 
UT 2009 Jan 28 and 30 visits, we may attempt to esti- 
mate the shape of the planet's transmission spectrum 
using the remaining transits. Although the evidence 
for spots in the final 8.0 fj,m transit on UT 2009 Feb 
2 is somewhat weaker, we choose to exclude it on the 
grounds that it displays some of the same behaviors (in- 
creased depth, larger than usual timing offset) as the 
more strongly affected 3.6 and 4.5 ^m transits imme- 
diately preceding it. If we then average the remain- 
ing three 8.0 //m depths, we find depths of [0.6694% ± 
0.0061%,0.6865% ± 0.0078%,0.6831% ± 0.0052%] at 3.6, 
4.5, and 8.0 fim, respectively. These three val ues are con- 
sistent with the near-IR transit depth from iPont et al.l 
(|2008bD of 0.6906% ± 0.0083% (1.1 - 1.9 ju m), as weU as 
the bes t-fit visible light transit depth from lBallard et all 
(l2010bD . 0.663% ± 0.014% (0.35 - 1.0 /xm). Ground- 
based data provide additional constraints in the near-IR, 
inclu ding a H ban d tran sit depth of 0.707% ± 0.019% 
from lAlonso et al.l 1)20081 ) and a Ks transit depth of 
0.64% ± 0.03% from ICaceres et all (l2009lf1 . both from 
individual transit observations. 

We fit these data us ing the retr i eval t echnique de- 
scribed in Madhusudh an fc Seagerl (j2009() . which ex- 
plores the parameter space of a one-dimensional, 
hydrogen-rich model atmosphere. We compute line- 
by-line radiative transfer with the assumption of hy- 
drostatic equilibrium and we use parametric prescrip- 
tions for the relative abundances of II2O, CH4, CO, 
CO2. We also include other dominant visible-light 

^ The best-fit planet-star radius ratio reported by these authors 
is inconsistent with their best-fit depth. We re-fit their data with 
an equivalent model and conclude that this discrepancy is most 
likely the result of a mistake in the reported value for the radius 
ratio, as our best-fit depth is a good match for the value stated in 
the paper. 



20 Knutson et al. 




-0.002 [ f t 1 

-0.08-0.06-0.04-0.02 0.00 0.02 0.04 0.06-0.06 -0.04 -0.02 0.00 0.02 0.04 0.06-0.08-0.06-0.04-0.02 0.00 0.02 0.04 0.06 

Time from Measured Transit Center [d] 



Figure 13. This plot shows the six transits observed in Jan/Feb. 2009. The left panel shows both 3.6 fj,m transits, the middle panel shows 
both 4.5 /im transits, and the right panel shows both 8.0 /im transits. The upper part of each panel overplots the normalized photometry 
for each visit (filled circles), with the first visit in black and the second visit in red, along with the best-fit transit light curves (solid lines) 
where the orbital inclination, a/i?*, and transit time have been allowed to vary freely for each individual transit. The lower panel takes the 
difference between the two light curves (black filled circles) and compares it to the difference between the best-fit transit solutions (solid 
black line). Note that even when all transit parameters are allowed to vary freely, it is not possible to reproduce the sharp features visible 
during ingress and egress in the lower left panel. 



and infrared opacity sources, including Na, K, H2-H2 
collision-induced absorption, a nd Rayleigh sca t tering . 
Our molecular li ne data are from lRothman et al.l ()2005D . 
iFreedman et al.l (j2p08), Fr eedman (pers onal commu- 
nication, 2009), iKarko schka fc Tomaskol (|2010f ). and 
Karkoschka (persona l co mmun ication, 20 11). The H2- 
H2 op acities are from lBorvsow et al.l (|1997| ) and lBorvsowl 
(|2002l ). We fix the pressure-temp erature jP-T) profile to 
the best-fit dayside profile from iStevenson et all ()2010D 
and iMadhusudhan &: Seageil ()201lD ; it is possible to ob- 
tain a marginally improved fit to these data if we allow 
the P-T profile to vary freely in the fit, but the differences 
are not significant. We find that the observations can be 
explained to within the 1-cr uncertainties by a methane- 
poor model (green line in Fig. I14p that contains mixing 
ratios of H2O = 1.0 x 10"^ CO = 1.0xlO-^ and CH4 = 
1.0 X 10~^; the data used in this fit appear to be incon- 
sistent with methane abundances > 10~^. This model 
also includes CO2 = 1.0 x 10~^, but the concentration 
of this molecule is less well constrained, as it is degen- 
erate with the CO abundance in the 4.5 /im band. We 
do not expect strong absorption due to atomic Na and 
K in this temperature regime ([Sharp fc Burrowa.2007,) . 



and we therefore adopt Na and K mixing ratios of 0.1 x 
solar abundances. If we compare the visible-light transit 
depth of 0.650% froni this model to the value reported by 
IBallard et all (|2010bn we find that it is consistent at the 
0.5g level. Model tra nsmission spectra for GJ 436b from 
iShabram et al.l ()2011[ ) , such as the rescaled model includ- 
ing higher-order hydrocarbons (model "g" in Shabram et 
al.) also provide a reasonably good match to these data. 

We can reduce the disagreement between the measured 
transit depths and the green model in the 1 — 2 /.tm wave- 
length range by introducing an opaque cloud layer at 50 
mbar (grey model in Fig. I14|) . However, such a cloud 
layer would be inconsi stent with the dayside emission 
spectrum measured bv IStevenson et all (|2010D unless it 
was optically thin in the center of the dayside hemi- 
sphere, or only intermittently present as discussed in 
M.l.ll We also note that occultations of spots and other 
features on the star will have a stronger effect on the 
measured transit depth at shorter wavelengths, and it 
is therefore possible that these measurements (several of 
which were derived from individual transit observations) 
are unreliable for our purposes here. 

Returning to the Spitzer data, we find that our con- 
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Figure 14. Comparison between measured transit depths (red circles) and model transmission spectra, where the transit depth is defined 
as the square of the best-fit planet- star radius ratio in eac h band. We i nclude p reviou sly published visible and near-IR t ransit depths from 
(in order of increasing wavelength) IBallard et aLl Il2010bl) . IPont et al.l Il2008bl1 . lAlonso'ct al. (2008), and ICaceres et al.l ll2009i) along with 
the Spitzer transit depths at 3.6, 4.5, and 8.0 /tm. Open circle s indicate Spitzer observations in which the planet appears to transit regions 
of non-uniform brightness on the star, as discussed in il4.1.3l Model transmission spectra include an atmosphere with reduced methane 
and enhanced CO abundances (Stevenson ct al. 2010) in green, a methane-rich model similar to that described in Bcaulicu et al. (20ll]) 
in blue, and a methane-poor model with an opaque cloud deck at pressures below 50 mbar in grey, which provides a better match to the 
visible-light transit depth from [Ballard et al.l (|20lOb V All models are calculated using t he methods d escribed in Madhusudh an fc Seageil 
((2009), where we fix the dayside P-T profiles to the nominal best-fit profile from Steven son et al.l (120 10 ) . Wc find that allowing the P-T 
profile to vary freely in our fits has a negligible effect on the agreement between the data and the green best-fit model. Colored green, blue, 
and grey circles indicate the predicted values for these models integrated over the bandpasses of the observations. 



elusions about the atmospheric composition are strongly 
dependent on our choice of which transit depths to in- 
clude in our analysis. We illustrate this with a blue 
model in Fig. [Ml which contains H2O and CH4 mix- 
ing ratios of 5.0 x 10^** each and no C O or CO2, and 
is coin pa rable to the model pr esented in iBeaulieu et al.l 
(|2011h . IBeaulieu etHI (|2011l ) excluded the shallower 
3.6 /xm transit on UT 2009 Jan 9 and kept the deeper 
3.6 fim UT 2009 Jan 28 and 8.0 fj.m UT 2009 Feb 2 
visits in their analysis, and as a result they concluded 
that the planet's transmission spectrum contained strong 
methane features, as illustrated by this blue model. They 
argue that the correction for the intrapixel effect is de- 
generate with the transit depth for the UT 2009 Jan 9 
visit and that this visit is therefore unreliable, but wc find 
that there is good overlap between the x and y positions 
spanned by the in-transit and out-of-transit data. We ob- 
tain transit depths that are consistent at the O.Ict level 
when we fit for our intrapixel sensitivity correction us- 
ing either the entire light curve or the out-of-transit data 
alone. Although our 3.6 and 8.0 /im transit depths are 



in good agreement with the values obtained by Beaulieu 
et al., our best-fit transit depth for the 4.5 fj,m UT 2009 
Jan 17 is 2.5a larger. We note that Beaulieu et al. al- 
low a/Ri, and b to vary individually for each transit, and 
that their values for these parameters from the Jan 17 
transit fit are outliers when compared to other visits; we 
conclude that this is likely the cause of their shallower 
best-fit radius ratio. Despite this disagreement, we find 
that if we include the same transits as Beaulieu et al. in 
our analysis, we also produce a transmission spectrum 
that is consistent with strong methane absorption. 

If, as we propose, occulted regions of non-uniform 
brightness on the surface of the star are responsible for 
the discrepancies in the 3.6 and 4.5 fim transit depths, it 
will be difficult to provide a definitive characterization of 
GJ 436b's transmission spectrum with broadband Spitzer 
photometry. Our analysis suggests that the atmosphere 
of GJ 436b is likely under-abundant in methane and 
ov er-abundant in CO, co nsiste nt with the conclusions 
of I Stevenson et all ()2010f ) and iMadhusudhan fc Seageii 
(|2011l) . but in order to reach these conclusions we have 



22 



Knutson et al. 



1.0006 




-0.05 0.00 
Time from Measured Eclipse Center [d] 



0.05 



Figure 15. Photometry for eleven 8 ^tm secondary eclipses (filled 
circles), with detector effects removed. Individual visits have been 
aligned using the best-fit transit ephemeris and assuming a con- 
stant offset for the secondary eclipse. The best-fit secondary eclipse 
light curve is overplotted (solid line), and residuals from this curve 
are shown in the lower panel. The period spanned by ingress and 
egress is denoted as a grey shaded region; we find no significant 
deviations from a model in which the planet has a uniform sur- 
face brightness. The dotted line shows the best-fit transit light 
curve, rescaled to match the depth of the secondary eclipse shown 
here. The longer ingress and egress for the transit are due to the 
increased planet-star distance and correspondingly higher impact 
parameter during this event, which occurs close to apastron. 

assumed that we have correctly identified and excluded 
all transits in which the planet occults active regions on 
the star. However, if the fractional spot coverage on the 
star is sufhciently high, it is possible that all transits are 
affected by these regions, in which case we cannot draw 
any robust conclusions about the shape of the planet's 
transmission spectrum. 

4.3. Dayside Emission Spectrum and Limits on 
Variability 

We can use the eleven secondary eclipse depths listed 
in Table [3] to study the properties of the planet's day- 
side atmosphere. We take the error-weighted aver- 
age of the eclipse depths and find a combined value 
of 0.0452% ± 0.0027%, consis tent with the value o f 
0.054% ± 0.008% reported by iStevenson et"all (f20Tot ). 
Next we construct a combined light curve incorporating 
all eleven secondary eclipse observations, shown in Fig. 
[TSl Fig. [in] shows the equivalent combined 8 /im tran- 
sit light curve for comparison. As a check wc fit these 
combined data with a secondary eclipse light curve and 
find that the best-fit eclipse depth agrees exactly with 
this error-weighted average from the individual eclipse 
fits. Because the strongest constraints on the relative 
abundances of methane and CO come from the 3.6 and 
4.5 /im eclipse measurements, we do not expect the re- 
duced 8 fim error bar to affect the conclusions reached by 
Stevenson et al. regarding these molecules. If we com- 
pare our results to the two models plotted in Fig. 2 of 
Stevenson et al., we find that the revised 8 /im eclipse 
depth is best-described by a cooler model with an effec- 
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Figure 16. Photometry for four 8 fj,m transits (filled circles), with 
detector effects removed and light curves aligned using the best-fit 
transit ephemeris. The best-fit transit light curve is overplotted 
(solid line), and residuals from this curve are shown in the lower 
panel. The period spanned by ingress and egress is denoted as a 
grey shaded region; we find no significant deviations from the ex- 
pected spherical planet model, as might be expected if the planet 
was significantly oblate. The dotted line shows the best-fit sec- 
ondary eclipse light curve, rescaled to match the depth of the tran- 
sit light curve where the limb-darkening is set to zero. The shorter 
ingress and egress for the secondary eclipse are due to the reduced 
planet-star distance and correspondingly lower impact parameter 
during this event, which occurs shortly before periastron. 

tive blackbody temperature of 790 K (defined as the tem- 
perature needed to match the total integrated flux at all 
wavelengths) and a modestly enhanced (30x higher) wa- 
ter abundance, rather than the hotter 860 K model with 
weaker water absorption. We also calculate a revised 
brightness temperature for the planet in the 8 fim band, 
defined as the temperature required to match the ob- 
served planet-star flux ratio in this bandpass assuming 
that the planet radiates as a blackbody. We use the pa- 
rameters in Table [5] and assume a Phoenix atmosphere 
model with an effective tem perature of 3585 K and \og{g) 
equal to 4.843 ()Torresl2007D for the star, and find that the 
planet has a best-fit brightness temperature of 740±16 K. 

Returning to Fig. I15[ we examine the residuals 
from our best-fit eclipse solution to search for evidence 
of deviations during ingress and egres s caused by a 
non- uniform day-side surfac e brightness (jWilliams et al.l 
[2006t iRauscher etldl [20071 ). The primary effect of a 
non-uniform brightne ss distribution is to shift the best- 
fit eclipse time (e.g.. [Agol et al.ll2010( ). but in this case 
uncertainties in estimates for GJ 436b's orbital eccen- 
tricity and longitude of periastron prevent us from de- 
tecting the small (< 1 minute) timing offsets expected 
from this effect. This timing offset will also display a 
small wavelength-dependence, due to variations in the 
brightness distribution as seen in different bandpasses, 
but this signal is likely to be too weak to detect by com- 
pari ng to the existing 3.6 and 5.8 /im eclipse observations 
from IStevenson et al.l (|2010D . Instead, we seek to deter- 
mine if the shape of the 8 fim eclipse ingress and egress 
can be used to constrain the planet's day-side brightness 
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Figure 17. The left panel shows an orbital diagram for the GJ 436 
system. Distances and radii are drawn to scale, and the location of 
periastron is marked by a dotted line. Grey shaded regions indicate 
the locations of the planet during the transit and secondary eclipse, 
where the viewer is assumed to be at the top of the plot. The right- 
hand panel shows five snapshots from a general circulation model 
for this planet as seen at different orbital phases by an observer on 
the Earth. 



distribution. We compare the eclipse light curves for a 
uniform surface brightness disk to that of a local equi- 
libirum model (i.e., one with the radiative time set to 
zero so that each reg ion of the planet is at its local equ i- 
librium temperature: lHansenll2008t iBurrows et al.l[2008l ). 
and find that the peak-to-trough residuals between these 
light curves is only 0.002%, if the eclipse depth is a free 
parameter. This is approximately a factor of ten smaller 
than our measurement errors, as demonstrated by the 
binned residuals in Fig. 1151 As we increase the amount 
of energy advected to the planet's night side using the 
models described in iCowan fc Agoll (|2011t ). the location 
of the hot spot on the planet's day side shifts away from 
the substellar point and the overall temperature contrast 
decreases. Because we are not sensitive to the timing off- 
set caused by the shifted hot spot, the only effect of this 
increased advection is to homogenize the planet's tem- 
peratures, producing light curves increasingly similar to 
the uniform disk light curves. 

4.3.1. A Variability Study for GJ 436b 

Tidal dissipation is expected to have driven GJ 436b 
into a pseudo-synchronous rotation state in which the 
planet's spin frequency is nearly commensurate with 
the planet's instantaneous orbital frequency at perias- 
tron. There are several competing theories of the pseu- 
dosyn chronization process (see, e.g. llvanov fc Pap alqizoul 
I2007f l. We adopt the expression given by Hut (1981D: 
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Figure 18. Predicted 8 ^m emission for GJ 436b (open circles) 
as a function of time, f rom g eneral circulation models described 
in lLangton fc LaughlinI 120081 ). The periodic modulation in flux is 
primarily due to the changing orbital geometry as we watch the 
hotter dayside rotate in and out of view, with a secondary effect 
caused by the heating and cooling of the atmosphere as the planet 
moves from periastron to apastron and back. The predicted fluxes 
during the secondary eclipse, as indicated by the black arrows, 
are nearly constant in time. The horizontal black line and grey 
shaded region indicate the average secondary eclipse depth and 
corresponding Icr uncertainty. 

For GJ 436b, this relation gives Pspin — 2.32 days, which 
yields a 19-day synodic period for the star as viewed from 
a fixed longitude on the planet. GJ 436b also experiences 
an 83% increase in incident flux during the 1.3-day in- 
terval between apoastron and periastron. 

We have computed simple hydrodynamical models 
to assess whether the asynchronous rotation and time- 
varying insolation are likely to generate atmospheric 
flows that are sufficiently chaotic to produce observ- 
able orbit-to-orbit variability in the secondary eclipse 
depths. Our two-dimensional hydrodynamical model 
contains three free parameters. The first, pa^m, is the at- 
mospheric pressure at the 8 ^m photosphere; the second, 
X, corresponds to the fraction of the incoming optical 
flux that is absorbed at or above the 8 fim photosphere; 
and the third, pb, corresponds to the pressure at the base 
of our modeled layer. We adopt parameter values of ps^m 
= 100 mbar, pi, = 4.0 bar and X = 1.0 for these models, 
which puts our models light curve in good agreement 
with GJ 436b's average 8 /xm secondary eclipse depth. 
The full details of the co mputational scheme are th e 
same as those adopted in iLangton fc LaughlinI (120081) , 
with updates as described in iLaughlin et al I (|2009D . A 
model photometric light curve is then obtained by inte- 
grating at each time step over the planetary hemisphere 
visible from Earth, where we assume that each patch of 
the planet radiates with a black-body spectrum corre- 
sponding to the local temperature. 

The model is run for a large number of orbits, and a 
quasi-steady state surface flow emerges. The tempera- 
ture structure of this flow as seen from an observer in 
the direction of Earth at five equally spaced intervals 
in the orbit is shown in Figure 1171 ^-nd the model light 
curve over these five orbits is shown in Figure 1181 Over 
the course of a single orbit, the 8 fiia planet-to-star fiux 
ratio varies nearly sinusoidally from AF/F = 0.033% to 
0.043%. The model's fiux at secondary eclipse agrees well 
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with the observed value, and varies by only 0.5% peak- 
to-peak from one orbit to the next. We note that more 
sophisticated three- dimensional genera l circulation mod- 
els for GJ 436b from lLewis et al.r(|201Cl[ ) also predict very 
low (1.3 — 1.5%) levels of variability in the 8 /im band 
for a range of atmospheric metallicities. 

Although these models indicate that GJ 436b's modest 
orbital eccentricity is likely not sufhcient to induce sig- 
nificant variability, they also do not include many pro- 
cesses such as clouds, photochemistry, and small scale 
turbulence that are known to contribute to temporal 
variability in planetary atmospheres. We therefore place 
empirical limits on GJ 436b's dayside variability using 
the eleven 8 fim secondary eclipse observations. We as- 
sume that the intrinsic dayside fluxes are drawn from 
either a Gaussian distribution with a standard deviation 
(5, or from a boxcar distribution with a width equal to 26. 
In both cases we set the mean of the distribution equal 
to the error-weighted mean of the measured secondary 
eclipse depths given in Table [5] We then conduct 10,000 
random trials, where we draw eleven measurements from 
each distribution and calculate the reduced of these 
values as compared to the measured secondary eclipse 
depths in Tabled We then determine the fraction of the 
10,000 random trials in which the reduced is less than 
or equal to one, which should correspond to the proba- 
bility that the underlying distribution is consistent with 
the measured eclipse depths. We repeat this calculation 
for a range of values for S, and plot the resulting proba- 
bility distribution as a function of 5 for both boxcar and 
Gaussian distributions. We find that for a boxcar distri- 
bution we can place [la, 2a, 3a] limits on the intrinsic 
variability of [29%, 42%, 58%], and for a Gaussian dis- 
tribution our corresponding upper limits are [17%, 27%, 
42%]. These limits are consistent with the predictions 
from general circulation models for this planet, but they 
are not low enough to provide meaningful constraints on 
these models. 

5. CONCLUSIONS 

In this paper we present Spitzer observations of eight 
transits and eleven secondary eclipses of GJ 436b at 3.6, 
4.5, and 8.0 /xm, which allow us to derive improved val- 
ues for the planet's orbital ephemeris, eccentricity, incli- 
nation, radius, and other system parameters. We discuss 
the effects that our assumptions about the longitude of 
periastron and stellar limb-darkening profiles have on our 
best-fit transit parameters, and find that our best-fit pa- 
rameters vary by la or less in all cases. We find that all 
parameters are consistent with a constant value over the 
two-year period spanned by our observations, with the 
exception of the measured transit depths and times in 
the 3.6 and 4.5 /im bands. We find that the 3.6 /im ra- 
dius ratio measured on UT 2009 Jan 28 is 4.7ct deeper 
than the value measured on UT 2009 Jan 9 in this same 
band, and the 4.5 ^m radius ratio from UT 2009 Jan 30 
is 2.9cr deeper than the value measured on UT 2009 Jan 
17. The level of significance for these changing radius 
ratios remains high even after accounting for the effects 
of residual correlated noise in the data. 

We also present an improved estimate for GJ 436b's 
8 /im secondary eclipse depth, based on eleven eclipse 
observations in this bandpass. We find that the new 
depth is consistent with previous models described 



in iStevenson et ahl ()2010[ ) and iMadhusudhan fc Seageii 
(|2011D . although we prefer solutions with modestly lower 
effective temperatures (790 K instead of 860 K). We use 
the shape of the eclipse ingress and egress to search for 
the presence of a non-uniform temperature distribution 
in the planet's dayside atmosphere, but uncertainties in 
the predicted time of secondary eclipse ultimately lim- 
its our ability to place meaningful constraints on this 
quantity. Our eclipse depths in this band are consistent 
with a constant value, and we place a Icr upper limit of 
17% on variability in the planet's dayside atmosphere. 
This limit is in good agreement with the predictions of 
general circulation models for this planet, which are typ- 
ically variable at the level of a few percent or less in this 
bandpass. 

Although it is possible that such residual noise or a 
time- varying cloud layer at sub-mbar pressures could ex- 
plain the apparent transit depth variations, the features 
observed in the transit light curves appear to be most 
consistent with the presence of occulted spots or other ar- 
eas of non-uniform brightness on the surface of the star in 
the UT 2009 Jan 28 and 30 transits. We find that for the 
UT 2009 Jan 28 transit the in-transit data have a higher 
RMS than the out-of-transit data, as would be expected 
for occulted spots; we would expect poorly corrected sys- 
tematics to produce an equivalently large RMS in both 
the in-transit and out-of-transit data, as the star spans 
same region of the pixel in both segments. Although we 
are not as sensitive to such effects in the UT 2009 Jan 30 
visit, which has higher levels of correlated noise due to an 
imperfect correction for intrapixel sensitivity variations, 
the short separation between these two observations rel- 
ative to the star's approximately 50 day rotation period 
means that the planet is likely to have occulted the same 
feature in both visits. We also see statistically significant 
variations in the measured transit times, where the am- 
plitude of the variations is typically smaller for infrared 
observations than for those obtained in visible light, also 
suggesting the present of occulted spots. We note that 
the anomalously deep transits observed on UT 2009 Jan 
28 and 30 also have best-fit transit times that are off- 
set by 30 s (3.1 — 3.5cr significance) from the predicted 
values. The fact that the three deepest transits are all 
measured within the same five-day period is also consis- 
tent with a single epoch of increased stellar activity. We 
reconcile this conclusion with the absence of any varia- 
tions larger than a few mmag in the star's visible and 
infrared fiuxes by proposing that the star's spin axis is 
likely inclined with respect to our line of sight, which has 
the effect of reducing the amplitude of any flux variations 
independent of spot coverage. If this is in fact the case, 
GJ 436b's orbit will be misaligned with respect to the 
star's spin axis. 

If we examine the wavelength-dependent transit depths 
for the subset of visits that appear to be least affected by 
spots, we find that the resulting transmission spectrum is 
consistent with the sam e reduced methane and enhanced 
CO abundances used bv IStevenson et aH (|2010f ) to fit the 
planet's dayside emission spectrum. These same tran- 
sit data are also consistent with models including an 
opaque cloud layer at a pressure of approximately 50 
mbar or less in the planet's atmosphere, which reduces 
the amplitude of the absorption features in the model 
spectra. We find no convincing evidence for the strong 
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methane absorption reported by iBeaulieu et all ()2011[ ). 
although we note that our conclusions vary significantly 
depending on which transits we include in our analy- 
sis. It is possible that all measured transit depths are 
alTected to varying degrees by stellar activity, in which 
case it may not be feasible to characterize the planet's 
transmission spectrum using broadband photometry ob- 
tained over multiple epochs. Because active regions oc- 
culted by the planet display a characteristic wavelength- 
dependence and also alter the local shape of the tran- 
sit light curve, high signal-to-noise grism spectroscopy of 
the transit over multiple epochs would help to resolve 
this issue. Such observations would also provide an in- 
dependent test of the reliability of the Spitzer transit 
data; if similar apparent depth variations were observed 
in other data sets, it would provide a strong argument 
against the hypothesis that the apparent depth variations 
in these data might be the result of poorly corrected in- 
strument effects. Lastly, grism spectroscopy could also 
be used to search for time-varying clouds at sub-mbar 
pressures, which should produce a featureless transmis- 
sion spectrum with a uniformly increased depth when 
present, as compared to the standard cloud- free trans- 
mission spectrum. 

As indicated by its rotation rate and Ca II H & K emis- 
sion, GJ 436 is an old and relatively quiet early M star. 
If the apparent transit depth variations we describe here 
are indeed due to the occultation of active regions on the 
star, as appears likely, we would expect similar features 
to occur frequently in the transit light curves of other 
planets orbiting M dwarfs at all activity levels. GJ 1214 
is currently the only other M star known to host a transit- 
ing planet, and has a similar 53-day rotation period and 
a modestly lower 3000 K effective temperature a s com- 
pared to GJ 436 (Charbonncau ct al. 2009; Bcrt a et al.l 
1201 Ih . A majority of the published data on this system 
are in the visible and near-infrared wavelengths where 
star spots should be prominent, and several recent pa- 
pers report the presence of occulted spots in a subset 
of transit observati ons (Be rta et al.ll20lTI : iCarter et al.l 
1201 H iKundurthv et al.ll201ll ) . Such spots might also ac- 
count for the apparent disagreement in measurements 
of the planet's infrared transmission spectrum, which 
some authors find to be featureless (|Bean et al.l 120101 : 
iDesert et al.ll2011bD. wh ile others detect absorption fea- 
tures ijCroll et al.ll20lil ). HD 189733b is currently the 
only other exoplanet with repeated Spitzer transit ob- 
servations in the same band; alt hough this planet or bits 
a relatively active K star fe.g.. iKnutson et al]|2010[ ). it 
exhibits much smaller variations in the meas ured transit 
depths an d times as com pared to GJ 436b (lAgol et al.l 
[2010; Des ert et al.ll2011al ). This is perhaps not surpris- 
ing, as the relative fractional spot coverage, spot sizes, 
and spot temperatures may well be qualitatively different 
on K stars and M stars. 
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